Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_sechiba/sechiba_hydrol_arch.f90 @ 7475

Last change on this file since 7475 was 2277, checked in by sebastiaan.luyssaert, 10 years ago
DEV: tested 200 years over Europe. Switched off the adaptation of allocation to water stress
File size: 139.6 KB

Line
1	!==================================================================================================================================
2	! MODULE : sechiba_hydrol_arch
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module calculates the amount of water that is available for the leaves for transpiration,
10	! accounting for effects of hydraulic architecture following (Hickler et al. 2006)
11	!!
12	!!\n DESCRIPTION : (1) Calculate soil water potential from soil water content as calculated in hydrol.f90
13	!! according to van Genuchten (1980). Soil water potential is weighted by the relative share of the
14	!! of theroots in each layer.
15	!! (2) Calculate pressure difference between leaves and soil (Whitehead 1998)
16	!! (3) Calculate root (Weatherly 1998), sapwood (Magnani et al. 2000, Tyree et al. 1994,
17	!! Sperry et al. 1998) and leave resistances (Sitch et al. 2003.
18	!! (4) Calculate the supply of water for transpiration by the leaves using Darcy's law (Whitehead 1998)
19	!!
20	!! RECENT CHANGE(S): None
21	!!
22	!! REFERENCE(S) : None
23	!!
24	!! SVN :
25	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_sechiba/sechiba.f90 $
26	!! $Date: 2013-05-06 14:14:36 +0200 (Mon, 06 May 2013) $
27	!! $Revision: 1295 $
28	!! \n
29	!_ ================================================================================================================================
30
31	MODULE sechiba_hydrol_arch
32
33	USE ioipsl
34
35	! modules used :
36	USE constantes
37	USE constantes_soil
38	USE pft_parameters
39	USE ioipsl_para
40	USE function_library, ONLY: wood_to_height_eff, biomass_to_lai, biomass_to_coupled_lai
41	USE qsat_moisture
42
43	IMPLICIT NONE
44
45	! Private and public routines
46
47	PRIVATE
48	PUBLIC hydrol_arch, jdiffuco_trans_co2
49
50
51	LOGICAL, SAVE :: l_first_sechiba = .TRUE. !! Flag controlling the intialisation (true/false)
52	!$OMP THREADPRIVATE(l_first_sechiba)
53
54	CONTAINS
55
56
57	!! ================================================================================================================================
58	!! SUBROUTINE : hydrol_arch
59	!!
60	!>\BRIEF Calculate the amount of water that is available for the leaves for transpiration,
61	! accounting for effects of hydraulic architecture following (Hickler et al. 2006)
62	!!
63	!!\n DESCRIPTION : (1) Calculate soil water potential from soil water content as calculated in hydrol.f90
64	!! according to van Genuchten (1980). Soil water potential is weighted by the relative share of the
65	!! roots in each layer.
66	!! (2) Calculate pressure difference between leaves and soil (Whitehead 1998)
67	!! (3) Calculate root (Weatherly 1998), sapwood (Magnani et al. 2000, Tyree et al. 1994, Sperry et al. 1998)
68	!! and leave resistances (Sitch et al. 2003).
69	!! (4) Calculate the supply of water for transpiration by the leaves using Darcy's law (Whitehead 1998)
70	!!
71	!! RECENT CHANGE(S): None
72	!!
73	!! MAIN OUTPUT VARIABLE(S): :: transpir_supply
74	!!
75	!! REFERENCE(S) : +++COMPLETE+++
76	!! Hickler et al. 2006
77	!! van Genuchten (1980)
78	!! Whitehead 1998
79	!! Weatherly 1998
80	!! Magnani et al. 2000
81	!! Tyree et al. 1994
82	!! Sperry et al. 1998
83	!! Sitch et al. 2003
84	!!
85	!! FLOWCHART :
86	!!
87	!_ ================================================================================================================================
88
89	SUBROUTINE hydrol_arch (kjpindex, circ_class_biomass, circ_class_n, temp_sol_new, &
90	temp_air, swc, transpir_supply, dtradia, njsc, biomass, ind, soiltile, &
91	vir_transpir_supply, veget)
92
93
94	!! 0. Variable declaration
95
96	!! 0.1 Input variables
97	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only (unitless)
98	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
99	INTEGER(i_std), DIMENSION(:), INTENT (in) :: njsc !! Index of the dominant soil textural class
100	!! in the grid cell used for hydrology
101	!! (1-nscm, unitless)
102	REAL(r_std),DIMENSION (:), INTENT (in) :: temp_sol_new !! New ground temperature (K)
103	REAL(r_std),DIMENSION (:), INTENT (in) :: temp_air !! Air temperature (K)
104	REAL(r_std),DIMENSION (:,:), INTENT (in) :: ind !! Number of individuals at the stand level
105	!! @tex $(m^{-2})$ @endtex
106	REAL(r_std),DIMENSION (:,:), INTENT (in) :: soiltile !! Fraction of each soiltile (-)
107	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: swc !! Volumetric soil water content
108	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: circ_class_n !! number of trees within a circumference class (tree m-2)
109	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass components of the model tree
110	!! within a circumference class
111	!! class @tex $(g C ind^{-1})$ @endtex
112	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! PFT total biomass
113	!! @tex $(gC m^{-2})$ @endtex
114	REAL(r_std),DIMENSION (:,:), INTENT (in) :: veget !! Fraction of the vegetation that covers the soil
115	!! also project canopy cover - calculated by Pgap
116
117
118	!! 0.2 Output variables
119
120	REAL(r_std), DIMENSION(:,:), INTENT (out) :: transpir_supply !! Supply of water for transpiration
121	!! @tex $(mm dt^{-1})$ @endtex
122	REAL(r_std), DIMENSION(:,:), INTENT(out) :: vir_transpir_supply !! Supply of water for transpiration
123	!! @tex $(mm dt^{-1})$ @endtex
124
125	!! 0.3 Modified variables
126
127
128	!! 0.4 Local variables
129
130	INTEGER(i_std) :: ipts,ivm,icirc !! indices
131	INTEGER(i_std) :: ibdl,istm !! indices
132	INTEGER(i_std) :: tmp_ncirc !! indices
133	REAL(r_std) :: R_root !! Hydraulic resistance of the roots
134	REAL(r_std) :: R_leaf !! Hydraulic resistance of the leaves
135	REAL(r_std) :: R_sap !! Hydraulic resistance of sapwood
136	REAL(r_std) :: R_shoot !! Hydraulic resistance of shoot corrected for temperature
137	REAL(r_std), DIMENSION(ncirc) :: height !! Vegetation height per circ_class_biomass (m)
138	REAL(r_std), DIMENSION(kjpindex,nbdl,nstm) :: phi_soil !! Soil water potential per soil layer per soiltile higher
139	!! than -5 MPa (threshold for hygroscopic water) (MPa)
140	REAL(r_std), DIMENSION(ncirc) :: circ_class_transpir_supply !! @tex $(m^{3}s^{-1}tree^{-1})$ @endtex
141	REAL(r_std), DIMENSION(ncirc) :: vir_circ_class_transpir_supply !! @tex $(m^{3}s^{-1}tree^{-1})$ @endtex
142	REAL(r_std) :: delta_P !! Pressure difference between leaves and soil
143	REAL(r_std) :: vir_delta_P !! Virtual pressure difference between leaves and soil
144	REAL(r_std), SAVE, DIMENSION(nbdl+1) :: z_soil !! Variable to store depth of the different
145	!! soil layers (m)
146	REAL(r_std), DIMENSION(nvm) :: k_sap_cav !! sapwood conductivity decreased by cavitation
147	REAL(r_std), DIMENSION(nbdl,nvm) :: root_dens !! root density per layer
148	REAL(r_std), DIMENSION(nvm) :: rpc !! scaling factor for integral
149	REAL(r_std),DIMENSION(kjpindex,nbdl,nstm) :: phi_soiltile !! soil water potentential per soil layer per soil tile
150	REAL(r_std),DIMENSION (kjpindex,nbdl,nstm) :: swc_adj !! Volumetric soil water content, adjusted when it is lower
151	!! than or equal to the residual soil water content
152	REAL(r_std),DIMENSION(kjpindex,nvm) :: phi_soilroot !! Sum of soil water potential in each soil layer weighted
153	!! by the share of the roots in each layer (MPa)
154	REAL(r_std),DIMENSION(kjpindex,nvm) :: vir_phi_soilroot !! Virtual transpir supply to be used in land cover change
155	!! Sum of soil water potential in each soil layer weighted
156	!! by the share of the roots in each layer (MPa)
157
158	!_ =================================================================================================================================
159
160	! Set Transpir_supply to a value that is definitely larger than transpir
161	! that way the plant will not experience water stress in the first
162	! day of the simulation when there Transpir_supply is calculated before
163	! there is any biomass available. Biomass is prescribed in
164	! stomate_prescribe.f90
165	! We need this line, if only for PFT type 1, which otherwise will
166	! not be initialized.
167	transpir_supply(:,:) = val_exp
168	phi_soiltile(:,:,:) = -un
169
170	!! 1. Calculate plant soil water potential
171
172
173	! Soil water potential for each soiltile in each layer from
174	! soil water content according to Van Genuchten (1980)
175
176	DO ipts = 1,kjpindex
177
178	!---TEMP---
179	IF (ld_hydrolarch) THEN
180	WRITE(numout,*) 'soil texture',njsc
181	!!$ WRITE(numout,*) 'veg_soiltext',veg_soiltex
182	WRITE(numout,*) 'swc', swc(ipts,:,:)
183	WRITE(numout,*)'soiltile',soiltile(ipts,:)
184	ENDIF
185	!----------
186
187	!! 2.1 Calculate soil water potential
188	! Swc calculated in hydrol.f90 cannot be lower than or equal
189	! to the residual soil water content for each soiltile.
190	DO ibdl = 1, nbdl
191
192	DO istm = 1,nstm
193
194	IF (swc(ipts,ibdl,istm) .LE. mcr_fao(istm)) THEN
195
196	!+++CHECK+++
197	!Explain why you do this
198	swc_adj(ipts,ibdl,istm) = mcr_fao(nstm) + .00001
199	!+++++++++++
200
201	ELSE
202
203	swc_adj(ipts,ibdl,istm) = swc(ipts,ibdl,istm)
204
205	IF (l_first_sechiba) THEN
206	swc_adj(ipts,ibdl,istm) = min_sechiba
207	ENDIF
208
209	ENDIF
210
211	ENDDO
212
213
214	!! 2. Initialize variables on first call
215	!!$ !---TEMP---
216	!!$ DO istm = 1,nstm
217	!!$ WRITE(numout,*) '2, ', ( swc_adj(ipts,ibdl,istm) - mcr_fao(njsc(ipts)) )
218	!!$ WRITE(numout,*) '3, ', ( mcs_fao(njsc(ipts)) - mcr_fao(njsc(ipts)) )
219	!!$ WRITE(numout,*) '4, ', ( swc_adj(ipts,ibdl,istm) - mcr_fao(njsc(ipts)) ) / &
220	!!$ ( mcs_fao(njsc(ipts)) - mcr_fao(njsc(ipts)) )
221	!!$ WRITE(numout,*) '5, ', ( -1 / (1-1/nvan_fao(njsc(ipts))) )
222	!!$ WRITE(numout,) 'line 1: ', - ( cte_grav rho_h2o * kilo_to_unit)
223	!!$ WRITE(numout,) 'line 2: ', ( 1/(avan_fao(njsc(ipts))kilo_to_unit) )
224	!!$ WRITE(numout,*) 'line 3: ', ( swc_adj(ipts,ibdl,istm) - mcr_fao(njsc(ipts)) ) / &
225	!!$ ( mcs_fao(njsc(ipts)) - mcr_fao(njsc(ipts)) )
226	!!$ WRITE(numout,*) 'line 4: ', ( -1 / (1-1/nvan_fao(njsc(ipts))) )
227	!!$ WRITE(numout,*) 'line 5: ', (1/nvan_fao(njsc(ipts)) )
228	!!$ WRITE(numout,*) 'swc: ', swc(ipts,ibdl,istm)
229	!!$ WRITE(numout,*) 'swc_adj: ', swc_adj(ipts,ibdl,istm)
230	!!$ ENDDO
231	!!$ !----------
232
233	! Soil water potential for each soiltile in each layer from
234	! soil water content according to Van Genuchten (1980). Convert
235	! unit avan from mm-1 to m-1. Convert units soil matric potential
236	! from m to MPa by multiplying with cte_grav and rho_h2o.
237	! Soiltile = fraction of each soil column
238	! (three soil columns=bare soil, tree, grass)
239	DO istm = 1,nstm
240
241	IF (swc_adj(ipts,ibdl,istm) .GT. min_stomate .AND. &
242	(swc_adj(ipts,ibdl,istm) - mcr_fao(njsc(ipts))) .GT. min_stomate) THEN
243
244	phi_soiltile(ipts,ibdl,istm) = &
245	- ( cte_grav * rho_h2o * kilo_to_unit * &
246	( 1/(avan_fao(njsc(ipts))kilo_to_unit) ) &
247	( ( (swc_adj(ipts,ibdl,istm) - mcr_fao(njsc(ipts)) ) / &
248	( mcs_fao(njsc(ipts)) - mcr_fao(njsc(ipts)) ) ) ** &
249	( -1 / (1-1/nvan_fao(njsc(ipts))) ) - 1) ** &
250	(1/nvan_fao(njsc(ipts)) ) ) / mega_to_unit
251
252	ELSE
253
254	phi_soiltile(ipts,ibdl,istm) = zero
255
256	ENDIF
257
258	ENDDO
259
260	! Threshold for soil water potential: it cannot be lower than the
261	! soil water potential for hygroscopic water (= -5 MPa, Larcher 2001)
262	phi_soil(ipts,ibdl,:) = MAX(phi_soiltile(ipts,ibdl,:),-5.)
263
264	! Calculate soil water potential per layer by multiplying with
265	! the fraction of each soiltile
266	!!$ phi_soil(ipts,ibdl) = SUM(soiltile(ipts,:) * phi_soiltile(ipts,ibdl,:))
267
268	ENDDO
269
270	! following initialisation of swc
271	l_first_sechiba = .FALSE.
272
273
274	!---TEMP---
275	IF (ld_hydrolarch) THEN
276	WRITE(numout,*) 'phi_soiltile', phi_soiltile(ipts,:,:)
277	WRITE(numout,*) 'phi_soil', phi_soil
278	ENDIF
279	!----------
280
281
282	!! 3. Calculate plant resistance
283
284	!! 3.1 Calculate root profile
285	! Weighting soil water potential with root profile
286	z_soil(1)= 0.
287	z_soil(2:nbdl+1) = diaglev(1:nbdl)
288
289	! Scaling factor for integration
290	rpc(:) = un / ( un - EXP( -z_soil(nbdl) / humcste(:) ) )
291
292	! Calculate root density per layer per pft
293	!+++CHECK+++
294	! what is the root density for layer one ?
295	DO ibdl = 1, nbdl ! Loop over # soil layers
296
297	root_dens(ibdl,:) = rpc(:) * &
298	( EXP( -z_soil(ibdl)/humcste(:)) - &
299	EXP( -z_soil(ibdl+1)/humcste(:)) )
300
301	ENDDO ! Loop over # soil layers
302	!++++++++++
303
304	!! 3.2 Calculate plant resistance to water transport
305	! Calculate plant resistance to water transport according to following
306	! the principle of an electric circuit. The code is largely based on
307	! Hickler et al. 2006
308	DO ivm = 2,nvm
309
310	!! 3.2.1 Calculate PFT-specific characteristics
311	IF (is_tree(ivm)) THEN
312
313	! Calculate for each circumference classes
314	tmp_ncirc = ncirc
315
316	! Calculate the effective height (m) from biomass. The effective
317	! value is used here because the water needs to be transported
318	! through both the belowground and aboveground biomass. Using the
319	! height (thus not the effective height) would only account for the
320	! height due to the aboveground biomass.
321	height(:) = wood_to_height_eff(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
322
323	! Calculate soil water potential weighted for root density
324	! Note that the second soiltile in swc is for trees.
325	!+++CHECK+++
326	! This implementation relies on two strong assumptions. First,
327	! roots are assumed to follow a rigid distribution irrespective
328	! of how the water is distributed in the soil. So even if on
329	! average lots of water is available in deeper layer, the modeled
330	! plants will grow most roots in the top layers. Second it is
331	! assumed that the root_potential is identical to the soil_potential
332	! With these assumptions it was found to be impossible to grow
333	! oak, beech and birch even in central Europe because those species
334	! experienced too much water stress. To overcome this problem
335	! a tuning factor has been introduced. In the future this tuning factor
336	! should be removed by a morec mechanistic model for root-soil
337	! interface.
338	phi_soilroot(ipts,ivm) = SUM(root_dens(:,ivm) * phi_soil(ipts,:,2)) + &
339	phi_soil_tune(ivm)
340	!+++++++++++
341
342	ELSE
343
344	! There are no circumference classes for grasses
345	! and crops. Biomasses are stored in first circumference
346	! class
347	tmp_ncirc = 1
348
349	! Calculate height (m) from biomass
350	height(1) = circ_class_biomass(ipts,ivm,1,ileaf,icarbon) * &
351	sla(ivm) * lai_to_height(ivm)
352
353	! Calculate soil water potential weighted for root density. Note
354	! that the third soiltile in swc is for grasses and crops.
355	!+++CHECK+++
356	! See comments above for trees
357	phi_soilroot(ipts,ivm) = SUM(root_dens(:,ivm) * phi_soil(ipts,:,3)) + &
358	phi_soil_tune(ivm)
359	!+++++++++++
360
361	ENDIF
362
363	!! 3.2.2 Calculate general characteristics
364	! Calculate soil water potential weighted for root density
365	! Note that the second soiltile in swc is for trees. Calculate
366	! for both the relevant soil water column i.e. tall vegetation
367	! and a vertual soil water colum i.e. bare soil. The latter is
368	! needed in land cover change because we may establish vegetation
369	! on PFT1 and then need to know the water stress the plants will
370	! experience on this new soil column.
371	!+++CHECK+++
372	! should this be sum of soil water potential or the average? The
373	! virtual_phi_soilroot may need to be revised once the water
374	! columns are getting merged in LCC. There are no plants so we
375	! should not account for the assumptions and thus not implement the
376	! tuning factor.
377	vir_phi_soilroot(ipts,ivm) = SUM(root_dens(:,ivm) * phi_soil(ipts,:,1))
378	!+++++++++++
379
380	!---TEMP---
381	IF (ld_hydrolarch .AND. ivm .EQ. test_pft) THEN
382	WRITE(numout,*) 'circ_class_biomass in hydrolic architecture,icirc', &
383	circ_class_biomass(1,test_pft,:,:,icarbon)
384	WRITE(numout,*) 'vegetation height, pft', height(:), ivm
385	WRITE(numout,*) 'z_soil',z_soil(1:nbdl+1)
386	WRITE(numout,*) 'rpc', ivm, rpc(ivm)
387	!!$ WRITE(numout,) 'phi_soilrootlayer',root_densphi_soil
388	WRITE(numout,*) 'phi_soilroot', phi_soilroot(:,test_pft)
389	WRITE(numout,*) 'vir_phi_soilroot', vir_phi_soilroot(:,test_pft)
390	ENDIF
391	!----------
392
393	!! 3.2.3 Calculate plant resistences
394	IF ( biomass(ipts,ivm,ileaf,icarbon) .GT. min_sechiba .AND. &
395	biomass(ipts,ivm,iroot,icarbon) .GT. min_sechiba ) THEN
396
397	DO icirc = 1,tmp_ncirc
398
399	! This ensures that the loop goes over all circ_classes for trees and only
400	! the first circ_class for grasses and crops
401	IF (circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon) .GT. min_sechiba .AND. &
402	circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) .GT. min_sechiba ) THEN
403
404	!! 3.2.3.1 Root resistance
405	! Root resistance is in the first order related to root mass
406	! and expression was proposed by Weatherly, 1982). Convert
407	! units of circ_class_biomss i.e. gC/ind to kg/ind
408	R_root = un / (k_root(ivm) * 2 * &
409	circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) &
410	/ kilo_to_unit)
411
412	!! 3.2.3.2 Sapwood resistance
413	! Calculate k_sap_cav as a function of water stress following
414	! Tyree et al. (1994) and Sperry et al. (1998)
415	IF (is_tree(ivm)) THEN
416
417	IF( (phi_soilroot(ipts,ivm) == zero) .AND. (phi_50(ivm) .LE. zero) )THEN
418
419	! This seems to cause a strange numerical issue, even though we should
420	! just be taking the exponential of zero (i.e., one).
421	k_sap_cav(ivm) = k_sap(ivm)
422
423	ELSE
424
425	k_sap_cav(ivm) = k_sap(ivm) * &
426	EXP(-((phi_soilroot(ipts,ivm)/phi_50(ivm)) ** &
427	c_cavitation(ivm)))
428	ENDIF
429
430	! If k_sap_cav is too low, we can get a divide by zero
431	! error in the next couple lines. So this is a protection
432	! against that.
433	IF(k_sap_cav(ivm) .LE. min_sechiba)THEN
434
435	k_sap_cav(ivm) = min_sechiba
436
437	ENDIF
438
439	! Following the expression by Magnani et al 2000, sapwood resistance
440	! is given as a function of sapwood area. Convert Sapwood mass to
441	! sapwood area
442	R_sap = height(icirc) / (k_sap_cav(ivm) * &
443	(( circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) + &
444	circ_class_biomass(ipts,ivm,icirc,isapbelow,icarbon)) / &
445	(tree_ff(ivm)*pipe_density(ivm) / height(icirc))) )
446
447	ELSE
448
449	!+++CHECK+++
450	! We have about the same leaf and root restistances and masses
451	! as for tress but lack a resistance in the sapwood. This is resulting
452	! in a huge supply because we lack part of the resistance. This was
453	! the motivation to create an R_sap for grasses and roots. This equation
454	! is made up and should be carefully checked against physiological
455	! studies for grasses and crops on the issue.
456	!!$ R_sap = zero
457	R_sap = height(icirc) / (k_sap(ivm) * &
458	(circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) / &
459	(pipe_density(ivm) / height(icirc))) )
460	!+++++++++++
461
462	ENDIF
463
464	!! 3.2.3.3 Leaf resistance
465	! Leaf resistance related to foliar projective area (Hickler 2006) or
466	! LAI (more consistent with pipe model)? fpc instead of lai because
467	! it is assumed that only sunlit leaves contribute to transpiration
468	! Initialy fpc = 1-e*(-0.5LAI) but it was replaxced by veget which
469	! in itseld is calculated in Pgap and thus accounts for LAI and canopy
470	! structure. This was still inconsistent with the calculation of the
471	! atmospheric demand for transpiration, to improve consistency the coupled
472	! lai was used as was done in the energy budget calculations.
473	! R_leaf = 1/(k_leaf * fpc)
474	!!$ R_leaf = un / (k_leaf(ivm) * &
475	!!$ (un - EXP(-0.5 * circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon) * &
476	!!$ sla(ivm))))
477	R_leaf = un / (k_leaf(ivm) * &
478	biomass_to_coupled_lai(circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon), &
479	veget(ipts,ivm),ivm))
480
481	!! 3.2.3.4 Temperature dependance of resistance
482	! Decrease root resistance with soil temperature as a result of
483	! changes in water viscosity (Cochard et al 2000)
484	! Decrease shoot resistance with air temperature (Cochard et al. 2000)
485	R_root = R_root / (a_viscosity(1) + a_viscosity(2) * &
486	(temp_sol_new(ipts) - zeroCelsius))
487	R_shoot = (R_leaf + R_sap) / (a_viscosity(1) + a_viscosity(2) * &
488	(temp_air(ipts) - zeroCelsius))
489
490	!! 3.2.3.5 Calculate pressure difference between leaves and soil
491	! Defined by Whitehead 1998. Convert units to MPa.
492	!+++CHECK+++
493	! If height is kept as a variable the delta_P gets out of bounds
494	! this seems to be one of the places where the static approach
495	! by Hickler fails. Could this be overcome by Sperry et al?
496	! As a quick fix we limit the height. This might have been acceptable
497	! if we would not have had to limit it at only 8m!
498	!!$ IF (height(icirc) .GT. 8.) height(icirc) = 8.
499	delta_P = phi_leaf(ivm) - phi_soilroot(ipts,ivm) + &
500	(height(icirc)cte_gravrho_h2o) / kilo_to_unit
501	vir_delta_P = phi_leaf(ivm) - vir_phi_soilroot(ipts,ivm) + &
502	(height(icirc)cte_gravrho_h2o) / kilo_to_unit
503	!++++++++++
504
505	!! 3.2.3.6 Calculate water supply for transpiration in m/s
506	! Using Dracy's law (Slatyer 1967, Whitehead 1998)
507	circ_class_transpir_supply(icirc) = (-delta_P) / (R_root + R_shoot)
508	vir_circ_class_transpir_supply(icirc) = (-vir_delta_P) / (R_root + R_shoot)
509
510	ELSE
511
512	!set transpir_supply to zero if there is no leaf or root biomass in the circ_class
513	circ_class_transpir_supply(icirc) = zero
514	vir_circ_class_transpir_supply(icirc) = zero
515
516	ENDIF
517
518	!---TEMP---
519	IF (ld_hydrolarch .AND. ivm .EQ. test_pft) THEN
520	WRITE(numout,*) 'root_dens', root_dens(:,test_pft)
521	WRITE(numout,*) 'k_sap_cavitation, biomass(iroot)', k_sap_cav(ivm), &
522	2 * circ_class_biomass(ipts,ivm,icirc,iroot,icarbon) / 1000.
523	WRITE(numout,*) 'sapwood area, ', ( circ_class_biomass(ipts,ivm,icirc,isapabove,icarbon) + &
524	circ_class_biomass(ipts,ivm,icirc,isapbelow,icarbon)) / (tree_ff(ivm) * pipe_density(ivm)) / &
525	height(icirc)
526	WRITE(numout,) 'fpc, ',(1-EXP(-0.5circ_class_biomass(ipts,ivm,icirc,ileaf,icarbon)*sla(ivm)))
527	WRITE(numout,*) 'R_root, R_sap, R_leaf', R_root, R_sap, R_leaf
528	WRITE(numout,*) 'R_shoot', R_shoot
529	WRITE(numout,*) 'delta_P', delta_P
530	ENDIF
531	!----------
532
533	ENDDO ! ncirc
534
535
536	!+++CHECK+++
537	! Just to be sure that we are not double counting for grasses and crops
538	IF (is_tree(ivm)) THEN
539	! transpir_supply per PFT (sum of diameter classes) per timestep
540	transpir_supply(ipts,ivm) = &
541	(SUM(circ_class_transpir_supply(:) * circ_class_n(ipts,ivm,:))) * &
542	dtradia * kilo_to_unit
543	vir_transpir_supply(ipts,ivm) = &
544	(SUM(vir_circ_class_transpir_supply(:) * circ_class_n(ipts,ivm,:))) * &
545	dtradia * kilo_to_unit
546	ELSE
547	transpir_supply(ipts,ivm) = &
548	(circ_class_transpir_supply(1) * circ_class_n(ipts,ivm,1)) * &
549	dtradia * kilo_to_unit
550	vir_transpir_supply(ipts,ivm) = &
551	(vir_circ_class_transpir_supply(1) * circ_class_n(ipts,ivm,1)) * &
552	dtradia * kilo_to_unit
553	ENDIF
554	!+++++++++++
555
556	! Set transpir_supply to zero when negative.
557	! (When phisoilroot becomes very low, delta_P becomes positive and transpir_supply negative)
558	IF (transpir_supply(ipts,ivm) .LT. zero) then
559	IF(ld_warn)THEN
560	WRITE(numout,*) 'transpir_supply before set to zero when negative', transpir_supply(ipts,ivm), ivm
561	WRITE(numout,*) 'ipts, ivm:', ipts,ivm
562	WRITE(numout,*) 'transpir_supply: ', transpir_supply(ipts,ivm)
563	WRITE(numout,*) 'phisoilroot: ', phi_soilroot(ipts,ivm)
564	ENDIF
565	transpir_supply(ipts,ivm) = MAX(transpir_supply(ipts,ivm), zero)
566	vir_transpir_supply(ipts,ivm) = MAX(vir_transpir_supply(ipts,ivm), zero)
567	ELSE
568	!do nothing
569	ENDIF
570
571	! transpir_supply is in m/s. Convert to mm/dt
572	! transpir_supply(ipts,ivm) = transpir_supply(ipts,ivm) * dtradia * 1000
573	!!$ WRITE(numout,*) 'transpir_supply', transpir_supply
574
575	!---TEMP---
576	IF (ld_hydrolarch .AND. ivm .EQ. test_pft .AND. ipts == test_grid) THEN
577	WRITE(numout,*) 'transpir_supply, PFT, ', transpir_supply(ipts,test_pft), ivm
578	WRITE(numout,*) 'vir_transpir_supply, PFT, ', vir_transpir_supply(ipts,test_pft), ivm
579	ENDIF
580	!----------
581
582	! Calculate daily water stress for use in stomate
583	!!$ IF (transpir_supply .LE. transpir) THEN
584	!!$ hydrol_stress(ipts,ivm) = transpir_supply(ipts,ivm)/transpir(ipts,ivm)
585	!!$ ELSE IF (transpir_supply .GT. transpir) THEN
586	!!$ hydrol_stress(ipts,ivm) = 1.
587	!!$ ENDIF
588
589	! Hydraulic architecure is not defined i.e. no biomass
590	ELSE
591
592	transpir_supply(ipts,ivm) = zero
593	phi_soilroot(ipts,ivm) = zero
594
595	ENDIF ! biomass leaves and roots.gt.zero
596
597
598	ENDDO ! nvm
599
600	ENDDO ! kjpindexipts
601
602	!---TEMP---
603	IF (ld_hydrolarch) THEN
604	WRITE(numout,*) 'transpir_supply, ', transpir_supply(:,:)
605	WRITE(numout,*) 'vir_transpir_supply, ', vir_transpir_supply(:,:)
606	ENDIF
607	!----------
608
609	END SUBROUTINE hydrol_arch
610
611
612
613	!! ==============================================================================================================================
614	!! SUBROUTINE : jdiffuco_trans_co2
615	!!
616	!>\BRIEF This subroutine computes carbon assimilation and stomatal
617	!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987).
618	!!
619	!! DESCRIPTION :\n
620	!! *** General:\n
621	!! The equations are different depending on the photosynthesis mode (C3 versus C4).\n
622	!! Assimilation and conductance are computed over 20 levels of LAI and then
623	!! integrated at the canopy level.\n
624	!! This routine also computes partial beta coefficient: transpiration for each
625	!! type of vegetation.\n
626	!! There is a main loop on the PFTs, then inner loops on the points where
627	!! assimilation has to be calculated.\n
628	!! This subroutine is called by diffuco_main only if photosynthesis is activated
629	!! for sechiba (flag STOMATE_OK_CO2=TRUE), otherwise diffuco_trans is called.\n
630	!! This subroutine is called at each sechiba time step by sechiba_main.\n
631	!! *** Details:
632	!! - Integration at the canopy level\n
633	!! \latexonly
634	!! \input{diffuco_trans_co2_1.1.tex}
635	!! \endlatexonly
636	!! - Light''s extinction \n
637	!! The available light follows a simple Beer extinction law.
638	!! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.\n
639	!! \latexonly
640	!! \input{diffuco_trans_co2_1.2.tex}
641	!! \endlatexonly
642	!! - Estimation of relative humidity of air (for calculation of the stomatal conductance)\n
643	!! \latexonly
644	!! \input{diffuco_trans_co2_1.3.tex}
645	!! \endlatexonly
646	!! - Calculation of the water limitation factor\n
647	!! \latexonly
648	!! \input{diffuco_trans_co2_2.1.tex}
649	!! \endlatexonly
650	!! - Calculation of temperature dependent parameters for C4 plants\n
651	!! \latexonly
652	!! \input{diffuco_trans_co2_2.2.tex}
653	!! \endlatexonly
654	!! - Calculation of temperature dependent parameters for C3 plants\n
655	!! \latexonly
656	!! \input{diffuco_trans_co2_2.3.tex}
657	!! \endlatexonly
658	!! - Vmax scaling\n
659	!! Vmax is scaled into the canopy due to reduction of nitrogen
660	!! (Johnson and Thornley,1984).\n
661	!! \latexonly
662	!! \input{diffuco_trans_co2_2.4.1.tex}
663	!! \endlatexonly
664	!! - Assimilation for C4 plants (Collatz et al., 1992)\n
665	!! \latexonly
666	!! \input{diffuco_trans_co2_2.4.2.tex}
667	!! \endlatexonly
668	!! - Assimilation for C3 plants (Farqhuar et al., 1980)\n
669	!! \latexonly
670	!! \input{diffuco_trans_co2_2.4.3.tex}
671	!! \endlatexonly
672	!! - Estimation of the stomatal conductance (Ball et al., 1987)\n
673	!! \latexonly
674	!! \input{diffuco_trans_co2_2.4.4.tex}
675	!! \endlatexonly
676	!!
677	!! RECENT CHANGE(S): N. de Noblet 2006/06
678	!! - addition of q2m and t2m as input parameters for the
679	!! calculation of Rveget
680	!! - introduction of vbeta23
681	!!
682	!! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular
683	!! concentration
684	!!
685	!! REFERENCE(S) :
686	!! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal
687	!! conductance and its contribution to the control of photosynthesis under
688	!! different environmental conditions, Prog. Photosynthesis, 4, 221â 224.
689	!! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis
690	!! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol.,
691	!! 19, 519â538.
692	!! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of
693	!! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78â90.
694	!! - Johnson, I. R., and J. Thornley (1984), A model of instantaneous and daily
695	!! canopy photosynthesis, J Theor. Biol., 107, 531 545
696	!! - McMurtrie, R.E., Rook, D.A. and Kelliher, F.M., 1990. Modelling the yield of Pinus radiata on a
697	!! site limited by water and nitrogen. For. Ecol. Manage., 30: 381-413
698	!! - Yin, X. and Struik, P. C., 2009. C3 and C4 photosynthesis models: An overview from the perspective
699	!! of crop modelling, NJAS - Wageningen Journal of Life Sciences, 57, 27--38.
700	!!
701	!! FLOWCHART : None
702	!! \n
703	!_ ================================================================================================================================
704
705	SUBROUTINE jdiffuco_trans_co2 (kjpindex, dtradia, swdown, temp_sol, pb, qsurf, q2m, t2m, &
706	temp_growth, rau, u, v, q_cdrag, vegstress, &
707	assim_param, Ca, &
708	veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, &
709	cimean, gsmean, gpp, vbeta23, Isotrop_Abs_Tot_p,&
710	Isotrop_Tran_Tot_p, lai_per_level, new_leaf_ci, &
711	hydrol_flag2, hydrol_flag3, gs_distribution, &
712	gs_diffuco_output, james_netcdf_flag, gstot_component, gstot_frac, warnings, veget, biomass)
713
714	!
715	!! 0. Variable and parameter declaration
716	!
717
718	!
719	!! 0.1 Input variables
720	!
721	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
722	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
723	REAL(r_std),DIMENSION (:), INTENT (in) :: swdown !! Downwelling short wave flux
724	!! @tex ($W m^{-2}$) @endtex
725	REAL(r_std),DIMENSION (:), INTENT (in) :: temp_sol !! Skin temperature (K)
726	REAL(r_std),DIMENSION (:), INTENT (in) :: pb !! Lowest level pressure (hPa)
727	REAL(r_std),DIMENSION (:), INTENT (in) :: qsurf !! Near surface specific humidity
728	!! @tex ($kg kg^{-1}$) @endtex
729	! N. de Noblet - 2006/06 - addition of q2m and t2m
730	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity
731	!! @tex ($kg kg^{-1}$) @endtex
732	! In off-line mode q2m and qair are the same.
733	! In off-line mode t2m and temp_air are the same.
734	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature (K)
735	! N. de Noblet - 2006/06 - addition of q2m and t2m
736	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
737	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! air density @tex ($kg m^{-3}$) @endtex
738	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
739	!! @tex ($m s^{-1}$) @endtex
740	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
741	!! @tex ($m s^{-1}$) @endtex
742	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
743	!! @tex ($m s^{-1}$) @endtex
744	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: assim_param !! min+max+opt temps (K), vcmax, vjmax for
745	!! photosynthesis
746	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
747	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: Ca !! CO2 concentration inside the canopy
748	!! @tex ($\mu mol mol^{-1}$) @endtex
749	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vegstress !! Soil moisture stress (0-1,unitless)
750	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum vegetation fraction of each PFT inside
751	!! the grid box (0-1, unitless)
752	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (-)
753
754	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index @tex ($m^2 m^{-2}$) @endtex
755	!! @tex ($m s^{-1}$) @endtex
756	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
757	!! @tex ($kg m^{-2}$) @endtex
758	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
759	!! @tex ($kg m^{-2}$) @endtex
760	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will
761	!! transpire (unitless)
762	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: Isotrop_Abs_Tot_p !! Absorbed radiation per level for photosynthesis
763	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: Isotrop_Tran_Tot_p !! Absorbed radiation per level for photosynthesis
764	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: lai_per_level !! This is the LAI per vertical level
765	!! @tex $(m^{2} m^{-2})$
766
767	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass @tex $(gC m^{-2})$ @endtex
768
769	LOGICAL, DIMENSION (kjpindex), INTENT (in) :: hydrol_flag2 !! flag that 'trips' the alternative energy budget for each grid square
770	LOGICAL, DIMENSION (kjpindex, nvm), INTENT (in) :: hydrol_flag3 !! flag that 'trips' the alternative energy budget for each grid square
771	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (in) :: gs_distribution
772	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (in) :: gs_diffuco_output
773	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (in) :: gstot_component
774	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (in) :: gstot_frac
775
776
777	LOGICAL, INTENT(in) :: james_netcdf_flag !! Flag that controls the netCDF output routine
778
779
780	!
781	!! 0.2 Output variables
782
783	!! 0.3 Modified variables
784	! JR these are all defined as inout variables as we skip over PFTs and grid squares for which the hydrol stress constraint
785	! does not apply. We can to use the values the values previously calculated (within diffuco) for these cases.
786	!
787	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration (unitless)
788
789	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3pot !! Beta for Potential Transpiration
790	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: rveget !! stomatal resistance of vegetation
791	!! @tex ($s m^{-1}$) @endtex
792	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex
793	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: cimean !! mean intercellular CO2 concentration
794	!! @tex ($\mu mol mol^{-1}$) @endtex
795	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: gsmean !! mean stomatal conductance to CO2 (umol m-2 s-1)
796	REAL(r_Std),DIMENSION (kjpindex,nvm), INTENT (inout) :: gpp !! Assimilation ((gC m^{-2} s^{-1}), total area)
797	REAL(r_std),DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) :: warnings !! A warning counter
798	!
799	!
800	REAL(r_std), DIMENSION (kjpindex,nvm,nlai) :: new_leaf_ci !! intercellular CO2 concentration (ppm)
801	!
802	!! 0.4 Local variables
803	!
804	REAL(r_std), DIMENSION(kjpindex, nvm) :: coupled_lai
805	REAL(r_std), DIMENSION(kjpindex, nvm) :: total_lai
806
807	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot) :: gs_output
808	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot) :: gstot_output
809	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: assimi_store
810
811	REAL(r_std), DIMENSION (kjpindex) :: jwind !! Wind module (m s^{-1})
812
813	REAL(r_std),DIMENSION (kjpindex,nvm) :: sum_rveget_rstruct
814
815	REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax !! maximum rate of carboxylation
816	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
817	INTEGER(i_std) :: ji, jv, ilev, limit_photo !! indices (unitless)
818	!REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest !! intercellular CO2 concentration at the lowest
819	! !! LAI level
820	! !! @tex ($\mu mol mol^{-1}$) @endtex
821	REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap !! relative water quantity in the water
822	!! interception reservoir (0-1,unitless)
823	REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex
824	! Assimilation
825	LOGICAL, DIMENSION(kjpindex) :: assimilate !! where assimilation is to be calculated
826	!! (unitless)
827	INTEGER(i_std) :: nia,inia,nina,inina,iainia !! counter/indices (unitless)
828	INTEGER(i_std), DIMENSION(kjpindex) :: index_assi,index_non_assi !! indices (unitless)
829	REAL(r_std), DIMENSION(kjpindex) :: vc2 !! rate of carboxylation (at a specific LAI level)
830	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
831	REAL(r_std), DIMENSION(kjpindex) :: vj2 !! rate of Rubisco regeneration (at a specific LAI
832	!! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex
833	REAL(r_std), DIMENSION(kjpindex) :: assimi !! assimilation (at a specific LAI level)
834	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
835	!! (temporary variables)
836	REAL(r_std), DIMENSION(kjpindex) :: gstop !! stomatal conductance to H2O at topmost level
837	!! @tex ($m s^{-1}$) @endtex
838	REAL(r_std), DIMENSION(kjpindex) :: gs !! stomatal conductance to CO2
839	!! @tex ($\mol m^{-2} s^{-1}$) @endtex
840
841
842	REAL(r_std), DIMENSION(kjpindex) :: gamma_star !! CO2 compensation point (ppm)
843	!! @tex ($\mu mol mol^{-1}$) @endtex
844
845	REAL(r_std), DIMENSION(kjpindex) :: air_relhum !! air relative humidity at 2m
846	!! @tex ($kg kg^{-1}$) @endtex
847	REAL(r_std), DIMENSION(kjpindex) :: VPD !! Vapor Pressure Deficit (kPa)
848	REAL(r_std), DIMENSION(kjpindex) :: water_lim !! water limitation factor (0-1,unitless)
849
850	REAL(r_std), DIMENSION(kjpindex) :: gstot !! total stomatal conductance to H2O
851	!! Final unit is
852	!! @tex ($m s^{-1}$) @endtex
853	REAL(r_std), DIMENSION(kjpindex) :: assimtot !! total assimilation
854	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
855	REAL(r_std), DIMENSION(kjpindex) :: Rdtot !! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1)
856	REAL(r_std), DIMENSION(kjpindex) :: gs_top !! leaf stomatal conductance to H2O across all top levels
857	!! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex
858	REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??)
859	!! @tex ($g g^{-1}$) @endtex
860	REAL(r_std), DIMENSION(kjpindex) :: T_Vcmax !! Temperature dependance of Vcmax (unitless)
861	REAL(r_std), DIMENSION(kjpindex) :: S_Vcmax_acclim_temp !! Entropy term for Vcmax
862	!! accounting for acclimation to temperature (J K-1 mol-1)
863	REAL(r_std), DIMENSION(kjpindex) :: T_Jmax !! Temperature dependance of Jmax
864	REAL(r_std), DIMENSION(kjpindex) :: S_Jmax_acclim_temp !! Entropy term for Jmax
865	!! accounting for acclimation to temperature (J K-1 mol-1)
866	REAL(r_std), DIMENSION(kjpindex) :: T_Rd !! Temperature dependance of Rd (unitless)
867	REAL(r_std), DIMENSION(kjpindex) :: T_Kmc !! Temperature dependance of KmC (unitless)
868	REAL(r_std), DIMENSION(kjpindex) :: T_KmO !! Temperature dependance of KmO (unitless)
869	REAL(r_std), DIMENSION(kjpindex) :: T_gamma_star !! Temperature dependance of gamma_star (unitless)
870	REAL(r_std), DIMENSION(kjpindex) :: vc !! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 mâ2 sâ1)
871	REAL(r_std), DIMENSION(kjpindex) :: vj !! Maximum rate of e- transport under saturated light (mumol CO2 mâ2 sâ1)
872	REAL(r_std), DIMENSION(kjpindex) :: Rd !! Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1)
873	REAL(r_std), DIMENSION(kjpindex) :: Kmc !! MichaelisâMenten constant of Rubisco for CO2 (mubar)
874	REAL(r_std), DIMENSION(kjpindex) :: KmO !! MichaelisâMenten constant of Rubisco for O2 (mubar)
875	REAL(r_std), DIMENSION(kjpindex) :: gb !! Boundary-layer conductance (mol mâ2 sâ1 barâ1)
876	REAL(r_std), DIMENSION(kjpindex) :: fvpd !! Factor for describing the effect of leaf-to-air vapour difference on gs (-)
877	REAL(r_std), DIMENSION(kjpindex) :: low_gamma_star !! Half of the reciprocal of Sc/o (bar bar-1)
878	REAL(r_std) :: N_Vcmax !! Nitrogen level dependance of Vcmacx and Jmax
879	REAL(r_std) :: fcyc !! Fraction of electrons at PSI that follow cyclic transport around PSI (-)
880	REAL(r_std) :: z !! A lumped parameter (see Yin et al. 2009) ( mol mol-1)
881	REAL(r_std) :: Rm !! Day respiration in the mesophyll (umol CO2 mâ2 sâ1)
882	REAL(r_std) :: Cs_star !! Cs -based CO2 compensation point in the absence of Rd (ubar)
883	REAL(r_std), DIMENSION(kjpindex) :: Iabs !! Photon flux density absorbed by leaf photosynthetic pigments (umol photon mâ2 sâ1)
884	REAL(r_std), DIMENSION(kjpindex) :: Jmax !! Maximum value of J under saturated light (umol eâ mâ2 sâ1)
885	REAL(r_std), DIMENSION(kjpindex) :: JJ !! Rate of eâ transport (umol eâ mâ2 sâ1)
886	REAL(r_std) :: J2 !! Rate of all eâ transport through PSII (umol eâ mâ2 sâ1)
887	REAL(r_std) :: VpJ2 !! eâ transport-limited PEP carboxylation rate (umol CO2 mâ2 sâ1)
888	REAL(r_std) :: A_1, A_3 !! Lowest First and third roots of the analytical solution for a general
889	!! cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
890	REAL(r_std) :: A_1_tmp, A_3_tmp !! Temporary First and third roots of the analytical solution for a general
891	!! cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
892	REAL(r_std) :: Obs !! Bundle-sheath oxygen partial pressure (ubar)
893	REAL(r_std) :: Cc !! Chloroplast CO2 partial pressure (ubar)
894	REAL(r_std) :: ci_star !! Ci -based CO2 compensation point in the absence of Rd (ubar)
895	REAL(r_std) :: a,b,c,d,m,f,j,g,h,i,l,p,q,r !! Variables used for solving the cubic equation (see Yin et al. (2009))
896	REAL(r_std) :: QQ,UU,PSI,x1,x2,x3 !! Variables used for solving the cubic equation (see Yin et al. (2009))
897
898	REAL(r_std),DIMENSION (kjpindex,nvm) :: lai_total !! total LAI for column in question
899
900	REAL(r_std) :: cresist !! coefficient for resistances (??)
901	REAL(r_std) :: lai_min !! mininum lai to calculate photosynthesis
902	!! coefficient for resistances (??)
903
904	LOGICAL,DIMENSION(kjpindex,nlevels_tot) :: top_level !! A flag to see if a given level is considered
905	!! a "top" level or not. Must be done this way
906	!! because there could be several "top" levels.
907	REAL(r_std) :: totlai_top !! The sum of the LAI in all the levels that we
908	!! declare as the "top".
909	INTEGER(i_std) :: lev_count !! How many levels we have already taken for
910	!! the "top"
911	INTEGER :: jij
912
913	! @defgroup Photosynthesis Photosynthesis
914	! @{
915	! 1. Preliminary calculations\n
916	REAL(r_std),DIMENSION (kjpindex,nvm) :: rveget_min !! Minimal Surface resistance of vegetation
917	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) &
918	:: light_tran_to_level !! Cumulative amount of light transmitted
919	!! to a given level
920	LOGICAL,DIMENSION (kjpindex) :: lskip !! A flag to indicate we've hit an unusual
921	!! case were we may have numerical instability.
922	INTEGER :: ipts
923	LOGICAL :: llight
924	LOGICAL :: llai
925	INTEGER,SAVE :: istep=0
926	!_ ================================================================================================================================
927
928
929	! JR 211013 define wind speed (this is defined in full in the diffuco module, but not here, so need to replicate it)
930	jwind(:) = SQRT (u(:)u(:) + v(:)v(:))
931	gs_output = 0.0d0
932	gstot_output = 0.0d0
933	lai_total = 0.0d0
934	assimi = 0.0d0
935	assimi_store = zero
936	!
937	DO ji=1,kjpindex
938	DO jv = 1, nvm
939	DO ilev=1,nlevels_tot
940	lai_total(ji,jv) = lai_total(ji,jv) + lai_per_level(ji,jv,ilev) !! total LAI for column in question
941	END DO ! ilev = 1, nlevels_tot
942	!IF (test_grid == 1 .AND. jv ==test_pft) THEN
943	! WRITE(numout, '(A20, ES20.4)') 'GET limited Vbeta3:', vbeta3(ji, jv)
944	!ENDIF
945	END DO ! jv = 1, nvm
946	END DO ! ji = 1, klpindex
947	!
948	! 1. Photosynthesis parameters
949	!
950	vcmax(:,:) = assim_param(:,:,ivcmax)
951	!
952	! @addtogroup Photosynthesis
953	! @{
954	!
955	! 1.2 Estimate relative humidity of air (for calculation of the stomatal conductance).\n
956	!! \latexonly
957	!! \input{jdiffuco_trans_co2_1.3.tex}
958	!! \endlatexonly
959	! @}
960	!
961	! N. de Noblet - 2006/06 - We use q2m/t2m instead of qair.
962	! CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
963	! air_relhum(:) = &
964	! ( qair(:) * pb(:) / (0.622+qair(:)*0.378) ) / &
965	! ( qsatt(:)pb(:) / (0.622+qsatt(:)0.378 ) )
966	! @codeinc
967
968	! WRITE(,) '181113a test, pre call to qsatcalc sechiba_hydrol_arch line 858; temp_sol is: ', temp_sol
969
970	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
971	air_relhum(:) = &
972	( qsurf(:) * pb(:) / (Tetens_1+qsurf(:)* Tetens_2) ) / &
973	( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) )
974
975
976	VPD(:) = ( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) ) &
977	- ( qsurf(:) * pb(:) / (Tetens_1+qsurf(:)* Tetens_2) )
978	! VPD is needed in kPa
979	! We limit the impact of VPD in the range of [0:6] kPa
980	! This seems to be the typical range of values
981	! that the vapor pressure difference can take, although
982	! there is a paper by Mott and Parkhurst, Plant, Cell & Environment,
983	! Vol. 14, pp. 509--519 (1991) which gives a value of 15 kPa.
984	VPD(:) = MAX(zero,MIN(VPD(:)/10.,6.))
985
986	! @endcodeinc
987	! N. de Noblet - 2006/06
988	!
989	!
990	! 2. beta coefficient for vegetation transpiration
991	!
992	rstruct(:,1) = rstruct_const(1)
993	rveget(:,:) = undef_sechiba
994	rveget_min(:,:) = undef_sechiba
995	!
996
997	! JR 051113 need to preserve the value of vbeta3 modified by supply term constraint
998	! in addition other values of vbeta3 are passed back from the earlier call to this
999	! routine in sechiba, so no need to initialise again. It is an (inout) variable in
1000	! this version of the routine.
1001	! vbeta3(:,:) = zero
1002
1003	vbeta3pot(:,:) = zero
1004	gsmean(:,:) = zero
1005	gpp(:,:) = zero
1006	light_tran_to_level(:,:,:)=un
1007	!
1008	cimean(:,1) = Ca(:)
1009	!
1010	! @addtogroup Photosynthesis
1011	! @{
1012	! 2. Loop over vegetation types\n
1013	! @}
1014	!
1015	DO jv = 2,nvm
1016	!
1017	! @addtogroup Photosynthesis
1018	! @{
1019	!
1020	! 2.1 Initializations\n
1021	!! \latexonly
1022	!! \input{diffuco_trans_co2_2.1.tex}
1023	!! \endlatexonly
1024	! @}
1025	!
1026	! beta coefficient for vegetation transpiration
1027	!
1028	rstruct(:,jv) = rstruct_const(jv)
1029	cimean(:,jv) = Ca(:)
1030	!
1031	!! mask that contains points where there is photosynthesis
1032	!! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points.
1033	!! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not
1034	!! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation,
1035	!! temperature and relative humidity).
1036	!! For the points where assimilation is not calculated, variables are initialized to specific values.
1037	!! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
1038	!! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes
1039	!! the nina points with no assimilation.
1040	!
1041	! set minimum lai value to test if photosynthesis should be calculated
1042	! 0.0001 is abitrary
1043	! could be externalised
1044	lai_min = 0.0001
1045	!
1046	nia=0
1047	nina=0
1048	index_assi=0
1049	index_non_assi=0
1050	!
1051	!
1052	DO ji=1,kjpindex
1053	!
1054	! This is tricky. We need a certain amount of LAI to have photosynthesis.
1055	! However, we divide our canopy up into certain levels, so even if we
1056	! have a total amount that is sufficient, there might not be enough
1057	! in any particular layer. So we check each layer to see if there
1058	! is one that has enough to give a non-zero photosynthesis. In
1059	! the worst-case scenario (small LAI in every other level), this becomes
1060	! our top level for the gs calculation. If gs_top_level is zero below,
1061	! we have a divide by zero issue.
1062	llai=.FALSE.
1063	DO ilev=1,nlevels_tot
1064	IF(lai_per_level(ji,jv,ilev) .GT. lai_min) llai=.TRUE.
1065	ENDDO
1066
1067	IF ( llai .AND. &
1068	( veget_max(ji,jv) .GT. min_sechiba ) ) THEN
1069
1070	! Here we need to do something different than in the previous model.
1071	! We explicitly need to check for absorbed light. This occurs
1072	! because our albedo in the two stream model is a function
1073	! of the solar angle. If the sun has set, we have no
1074	! light, and therefore we cannot have photosynthesis.
1075	llight=.FALSE.
1076	DO ilev=1,nlevels_tot
1077	IF(Isotrop_Abs_Tot_p(ji,jv,ilev) .GT. min_sechiba) llight=.TRUE.
1078	ENDDO
1079	IF ( ( swdown(ji) .GT. min_sechiba ) .AND. &
1080	( temp_growth(ji) .GT. tphoto_min(jv) ) .AND. &
1081	llight .AND. &
1082	( temp_growth(ji) .LT. tphoto_max(jv) )) THEN
1083	!
1084	assimilate(ji) = .TRUE.
1085	nia=nia+1
1086	index_assi(nia)=ji
1087	!
1088	ELSE
1089	!
1090	assimilate(ji) = .FALSE.
1091	nina=nina+1
1092	index_non_assi(nina)=ji
1093	!
1094	ENDIF
1095	ELSE
1096	!
1097	assimilate(ji) = .FALSE.
1098	nina=nina+1
1099	index_non_assi(nina)=ji
1100	!
1101	ENDIF
1102	!
1103	ENDDO
1104	!
1105
1106	IF(ld_photo .AND. jv == test_pft)THEN
1107	WRITE(numout,*) 'index_non_assi(:): ',index_non_assi(:)
1108	WRITE(numout,*) 'index_assi(:): ',index_assi(:)
1109	WRITE(numout,*) 'nia: ',nia
1110	ENDIF
1111
1112
1113	top_level(:,:)=.FALSE.
1114	gstot(:) = zero
1115	assimtot(:) = zero
1116	Rdtot(:)=zero
1117	gs_top(:) = zero
1118	lskip(:)=.FALSE.
1119	!
1120	zqsvegrap(:) = zero
1121	WHERE (qsintmax(:,jv) .GT. min_sechiba)
1122	!! relative water quantity in the water interception reservoir
1123	zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
1124	ENDWHERE
1125	!
1126	!! Calculates the water limitation factor.
1127	!!Not used when hydrol_arch is used, can be deleted
1128	WHERE ( assimilate(:) )
1129	! water_lim(:) = MAX( min_sechiba, MIN( vegstress(:,jv)/0.5, un ))
1130	water_lim(:) = un
1131	ELSEWHERE
1132	water_lim(:) = zero
1133	ENDWHERE
1134
1135	! give a default value of ci for all pixel that do not assimilate
1136	DO ilev=1,nlevels_tot
1137	DO inina=1,nina
1138	new_leaf_ci(index_non_assi(inina),jv,ilev) = Ca(index_non_assi(inina))
1139	ENDDO
1140	ENDDO
1141	!
1142	!
1143	! Here is the calculation of assimilation and stomatal conductance
1144	! based on the work of Farquahr, von Caemmerer and Berry (FvCB model)
1145	! as described in Yin et al. 2009
1146	! Yin et al. developed a extended version of the FvCB model for C4 plants
1147	! and proposed an analytical solution for both photosynthesis pathways (C3 and C4)
1148	! Photosynthetic parameters used are those reported in Yin et al.
1149	! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use
1150	! Medlyn et al. (2002) and Kattge & Knorr (2007)
1151	! Because these 2 references do not consider mesophyll conductance, we neglect this term
1152	! in the formulations developed by Yin et al.
1153	! Consequently, gm (the mesophyll conductance) tends to the infinite
1154	! This is of importance because as stated by Kattge & Knorr and Medlyn et al.,
1155	! values of Vcmax and Jmax derived with different model parametrizations are not
1156	! directly comparable and the published values of Vcmax and Jmax had to be standardized
1157	! to one consistent formulation and parametrization
1158
1159	! See eq. 6 of Yin et al. (2009)
1160	! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
1161	T_KmC(:) = Arrhenius(kjpindex,temp_sol,298.,E_KmC(jv))
1162	T_KmO(:) = Arrhenius(kjpindex,temp_sol,298.,E_KmO(jv))
1163	T_gamma_star(:) = Arrhenius(kjpindex,temp_sol,298.,E_gamma_star(jv))
1164
1165
1166	! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
1167	T_Rd(:) = Arrhenius(kjpindex,temp_sol,298.,E_Rd(jv))
1168
1169
1170	! For C3 plants, we assume that the Entropy term for Vcmax and Jmax
1171	! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10
1172	! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function
1173	! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002)
1174	! and Kattge & Knorr (2007).
1175	! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function
1176	! Concerning this apparent unconsistency, have a look to the section 'Limitation of
1177	! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation
1178
1179	! Growth temperature tested by Kattge & Knorr range from 11 to 35Â°C
1180	! So, we limit the relationship between these lower and upper limits
1181	S_Jmax_acclim_temp(:) = aSJ(jv) + bSJ(jv) * MAX(11., MIN(temp_growth(:),35.))
1182	T_Jmax(:) = Arrhenius_modified(kjpindex,temp_sol,298.,E_Jmax(jv),D_Jmax(jv),S_Jmax_acclim_temp)
1183
1184	S_Vcmax_acclim_temp(:) = aSV(jv) + bSV(jv) * MAX(11., MIN(temp_growth(:),35.))
1185	T_Vcmax(:) = Arrhenius_modified(kjpindex,temp_sol,298.,E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
1186	!
1187	vc(:) = vcmax(:,jv) * T_Vcmax(:)
1188	! As shown by Kattge & Knorr (2007), we make use
1189	! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants
1190	! rJV is written as a function of the growth temperature
1191	! rJV = arJV + brJV * T_month
1192	! See eq. 10 of Kattge & Knorr (2007)
1193	! and Table 3 for Values of arJV anf brJV
1194	! Growth temperature is monthly temperature (expressed in Â°C) - See first paragraph of
1195	! section Methods/Data of Kattge & Knorr
1196	vj(:) = ( arJV(jv) + brJV(jv) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:,jv) * T_Jmax(:)
1197	!
1198	! @endcodeinc
1199	!
1200	KmC(:)=KmC25(jv)*T_KmC(:)
1201	KmO(:)=KmO25(jv)*T_KmO(:)
1202	gamma_star(:) = gamma_star25(jv)*T_gamma_star(:)
1203	! low_gamma_star is defined by Yin et al. (2009)
1204	! as the half of the reciprocal of Sco - See Table 2
1205	! We derive its value from Equation 7 of Yin et al. (2009)
1206	! assuming a constant O2 concentration (Oi)
1207	low_gamma_star(:) = gamma_star(:)/Oi
1208	! VPD expressed in kPa
1209	fvpd(:) = 1. / ( 1. / (a1(jv) - b1(jv) * VPD(:)) - 1. )
1210	! * water_lim(:)
1211	!???
1212	!leaf boundary layer conductance
1213	gb(:) = (1.0_r_std /25.0_r_std)/(22.4_r_std *temp_sol(:)/273._r_std/1000._r_std)
1214	!???
1215	!
1216	! @addtogroup Photosynthesis
1217	! @{
1218	!
1219	! 2.4 Loop over LAI levels to estimate assimilation and conductance\n
1220	! @}
1221	!
1222
1223
1224
1225
1226	IF (nia .GT. 0) THEN
1227	gridloop1: DO inia=1,nia
1228	!
1229	iainia=index_assi(inia)
1230
1231	! JR we 'cycle' the LAI loop for those PFTs not subject to the hydrol stress for the grid square in question
1232
1233	IF (hydrol_flag3(iainia, jv)) THEN
1234	!nothing
1235	ELSE
1236	CYCLE gridloop1
1237	END IF !(hydrol_flag3)
1238
1239
1240	! JR the first part of the routine is a copy of diffuco_trans_co2, so that we can
1241	! calculate values for fvpd, assimtot, Rdtot, new_leaf_ci and ci_star
1242	!
1243	! Find the top layer that still contains lai. We have also
1244	! added a check in here to make sure this layer has absorbed
1245	! light. If there is no absorbed light, we will not calculate
1246	! stomatal conductance for the top level, which means gstop
1247	! will be zero later and we'll crash when we divide by it.
1248	! Include a warning so we know how often this happens.
1249	! Include a warning so we know how often this happens.
1250	! The old system of chosing just one top level caused some bad results
1251	! in the transpir, evidently due to the calculation of vbeta3 being
1252	! too low due to too low of LAI in the top level. Therefore, we are
1253	! now going to take several top levels in an effort to make sure
1254	! we have enough LAI. We aim for an arbitrary amount of 0.1 or
1255	! the top three levels. This has not been rigoursly tested, but it
1256	! is theoretically justified by the fact that if the LAI of the top
1257	! level is very low, other leaves will be getting a lot of sun and
1258	! well-exposed to the atmosphere, and therefore contribute equal
1259	! to the resistence of the vegetation.
1260	totlai_top=zero
1261	lev_count=0
1262	!
1263	DO ilev=nlevels_tot, 1, -1
1264	IF(lai_per_level(iainia,jv,ilev) .GT. lai_min)THEN
1265	IF(Isotrop_Abs_Tot_p(iainia,jv,ilev) .GT. min_sechiba)THEN
1266	IF(totlai_top .LT. lai_top(jv) .AND. lev_count .LT. nlev_top)THEN
1267	lev_count=lev_count+1
1268	top_level(iainia,ilev)=.TRUE.
1269	totlai_top=totlai_top+lai_per_level(iainia,jv,ilev)
1270	ELSE
1271	! We have enough levels now.
1272	!WRITE(numout,*) 'Top levels ',totlai_top,lev_count
1273	EXIT
1274	ENDIF
1275	ELSE
1276	IF(ld_warn)THEN
1277	WRITE(numout,*) 'WARNING: Found a level with LAI, ' // &
1278	'but no absorbed light.'
1279	WRITE(numout,*) 'WARNING: iainia,jv,ilev, ',iainia,jv,ilev
1280	WRITE(numout,*) 'WARNING: Isotrop_Abs_Tot_p(iainia,jv,ilev),' //&
1281	'lai_per_level(iainia,jv,ilev), ',iainia,jv,ilev
1282	ENDIF
1283	ENDIF
1284	ENDIF
1285	ENDDO
1286	IF(totlai_top .LT. min_sechiba)THEN
1287	! There seem to be some fringe cases where this happens, usually
1288	! at sunrise or sunset with low LAI levels. Instead of stopping
1289	! here, I'm going to issue a warning. If this case happens,
1290	! we produce no GPP, but we're not going to crash. We need to
1291	! be careful in one of the loops below, since a stomatal
1292	! conductance of zero will cause a crash.
1293	IF(ld_warn)THEN
1294	WRITE(numout,*) 'WARNING: Did not find a level with LAI in diffuco!'
1295	WRITE(numout,*) 'WARNING: iainia,ivm: ',iainia,jv
1296	ENDIF
1297	assimtot(iainia) = zero
1298	Rdtot(iainia) = zero
1299	gstot(iainia) = zero
1300	lskip(iainia)=.TRUE.
1301	CYCLE gridloop1
1302	ENDIF
1303
1304	! We need to know how much is transmitted to this level
1305	! from the top of the canopy, since this aboslute number
1306	! is what is actually used in the photosyntheis. Right now
1307	! Isotrop_Tran_Tot_p and Isotropic_Abs_Tot_p are the
1308	! relative quantities transmitted through and absorbed
1309	! by the layer.
1310	DO ilev = nlevels_tot-1,1,-1
1311	light_tran_to_level(iainia,jv,ilev)=light_tran_to_level(iainia,jv,ilev+1)*Isotrop_Tran_Tot_p(iainia,jv,ilev+1)
1312	ENDDO
1313	IF(ld_photo .AND. jv == test_pft .AND. iainia == test_grid)THEN
1314	WRITE(numout,*) 'Absolute transmitted light to each level',jv,iainia
1315	DO ilev = nlevels_tot,1,-1
1316	WRITE(numout,*) ilev,light_tran_to_level(iainia,jv,ilev)
1317	ENDDO
1318	ENDIF
1319
1320	DO ilev = 1, nlevels_tot
1321
1322	IF (Isotrop_Abs_Tot_p(iainia,jv,ilev) .GT. min_sechiba) THEN
1323
1324	IF(lai_per_level(iainia,jv,ilev) .LE. min_sechiba)THEN
1325	! This is not necessarily a problem. LAI is updated in slowproc, while
1326	! absorbed light is updated in condveg, which is called after diffuco.
1327	! So we might be one timestep off. It generally should be at night,
1328	! though, so there is no absorbed light. If this persists, that
1329	! is a real problem.
1330	IF(ld_warn)THEN
1331	WRITE(numout,*) 'WARNING: We have absorbed light but no LAI ' // &
1332	'in this level. Skipping.'
1333	WRITE(numout,'(A,3I8)') 'WARNING: iainia,jv,ilev: ',iainia,jv,ilev
1334	WRITE(numout,'(A,2F20.10)') 'WARNING: Isotrop_Abs_Tot_p(iainia,jv,ilev), ' // &
1335	'lai_per_level(ji,jv,ilev): ',&
1336	Isotrop_Abs_Tot_p(iainia,jv,ilev), lai_per_level(iainia,jv,ilev)
1337	ENDIF
1338	CYCLE
1339	ENDIF
1340	!
1341	! nitrogen is scaled according to the transmitted light
1342	! before the caclulation was (1 - exp(lai)) we replaced
1343	! it by the transmitted light calculated by the 2stream
1344	! model, maybe ( un - Isotrop_Tran_Tot_p(iainia,jv,ilev) )
1345	! could be replaced by the directly calculated absorbed
1346	! light Isotrop_Abs_Tot_p(iainia,jv,ilev).
1347	! Notice that light in the old code was the cumulative light transmitted to
1348	! the layer, while Isotrop_Tran_Tot_p is the relative light transmitted
1349	! through the layer. Therefore, we use a different variable.
1350	!
1351	! N_Vcmax = ( un - .7_r_std * ( un - Isotrop_Tran_Tot_p(iainia,jv,ilev) ) )
1352	N_Vcmax = ( un - .7_r_std * ( un - light_tran_to_level(iainia,jv,ilev) ) )
1353	!
1354	! IF (control%ok_hydrol_arch) THEN
1355
1356	!*******CHECK: only photsynthesis when day*************
1357	! WHERE (assimilate(:))
1358	! vc2(iainia) = vc(iainia) * N_Vcmax
1359	! vj2(iainia) = vj(iainia) * N_Vcmax
1360	! see Comment in legend of Fig. 6 of Yin et al. (2009)
1361	! Rd25 is assumed to equal 0.01 Vcmax25
1362	! Rd(iainia) = vcmax(iainia,jv) * N_Vcmax * 0.01 * T_Rd(iainia)
1363	!ELSEWHERE
1364	! vc2(iainia) = zero
1365	! vj2(iainia) = zero
1366	! Rd(iainia) = zero
1367	!ENDWHERE
1368
1369	! ELSE
1370
1371	vc2(iainia) = vc(iainia) * N_Vcmax
1372	vj2(iainia) = vj(iainia) * N_Vcmax
1373
1374	! see Comment in legend of Fig. 6 of Yin et al. (2009)
1375	! Rd25 is assumed to equal 0.01 Vcmax25
1376	Rd(iainia) = vcmax(iainia,jv) * N_Vcmax * 0.01 * T_Rd(iainia)
1377
1378	!ENDIF
1379
1380	!+++++ CHECK +++++
1381	IF(jv == test_pft .AND. iainia == test_grid .AND. ld_photo)THEN
1382	WRITE(numout,*) 'Iabs1: ',iainia,jv,ilev
1383	WRITE(numout,'(A,10F20.10)') 'Iabs2: ',swdown(iainia),W_to_mmol,RG_to_PAR,&
1384	Isotrop_Abs_Tot_p(iainia,jv,ilev),light_tran_to_level(iainia,jv,ilev) &
1385	, lai_per_level(iainia,jv,ilev)
1386	ENDIF
1387	Iabs(iainia)=swdown(iainia)W_to_mmolRG_to_PAR*&
1388	Isotrop_Abs_Tot_p(iainia,jv,ilev)*light_tran_to_level(iainia,jv,ilev) &
1389	/ lai_per_level(iainia,jv,ilev) ! need to convert to per m**2 of leaves to
1390	! be consistent with the units we use for Jmax
1391	! WRITE(numout,*) 'Iabs ',jv,Iabs(iainia),lai_per_level(iainia,jv,ilev)
1392	! Iabs(iainia)=swdown(iainia)W_to_mmolRG_to_PAR*Isotrop_Abs_Tot_p(iainia,jv,ilev)
1393	!+++++++++++++++++
1394	! If we have a very dense canopy, the lower levels might
1395	! get almost no light. If the light is close to zero, it
1396	! causes a numerical problem below, as x1 will be zero. Of
1397	! course, if the absorbed light is really small photosynthesis will
1398	! not happen, so we can just skip this level. This number is
1399	! fairly arbitrary, but it should be as small as possible.
1400	IF(Iabs(iainia) .LT. 1e-12)CYCLE
1401
1402	!
1403	!
1404	! eq. 4 of Yin et al (2009)
1405	Jmax(iainia)=vj2(iainia)
1406	JJ(iainia) = ( alpha_LL(jv) * Iabs(iainia) + Jmax(iainia) - &
1407	sqrt((alpha_LL(jv) * Iabs(iainia) + Jmax(iainia) )**2. &
1408	- 4 * theta(jv) * Jmax(iainia) * alpha_LL(jv) * &
1409	Iabs(iainia)) ) &
1410	/ ( 2 * theta(jv))
1411
1412	!
1413
1414	IF(ld_photo .AND. test_pft == jv)THEN
1415	WRITE(numout,) '************************'
1416	WRITE(numout,*) 'jv,ilev: ',jv,ilev
1417	WRITE(numout,*) 'JJ(iainia),Jmax(iainia),Iabs(iainia),Rd(iainia): ',&
1418	JJ(iainia),Jmax(iainia),Iabs(iainia),Rd(iainia)
1419	WRITE(numout,*) 'vj2(iainia),vc2(iainia),N_Vcmax:',&
1420	vj2(iainia),vc2(iainia),N_Vcmax
1421	WRITE(numout,*) 'Isotrop_Abs_Tot_p(iainia,jv,ilev),Isotrop_Tran_Tot_p(iainia,jv,ilev):',&
1422	Isotrop_Abs_Tot_p(iainia,jv,ilev),Isotrop_Tran_Tot_p(iainia,jv,ilev)
1423
1424	WRITE(numout,*) 'water_lim(iainia),T_Rd(iainia):',water_lim(iainia),T_Rd(iainia)
1425	WRITE(numout,*) 'vc(iainia),vj(iainia), vcmax(iainia,jv):',&
1426	vc(iainia),vj(iainia),vcmax(iainia,jv)
1427	WRITE(numout,) '************************'
1428	ENDIF
1429	IF ( is_c4(jv) ) THEN
1430	!
1431	! @addtogroup Photosynthesis
1432	! @{
1433	!
1434	! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n
1435	!! \latexonly
1436	!! \input{jdiffuco_trans_co2_2.4.2.tex}
1437	!! \endlatexonly
1438	! @}
1439	!
1440	!
1441	!
1442	! Analytical resolution of the Assimilation based Yin et al. (2009)
1443	! Eq. 28 of Yin et al. (2009)
1444	fcyc= 1. - ( 4.(1.-fpsir(jv))(1.+fQ(jv)) + 3.h_protons(jv)fpseudo(jv) ) / &
1445	( 3.h_protons(jv) - 4.(1.-fpsir(jv)))
1446
1447	! See paragraph after eq. (20b) of Yin et al.
1448	Rm=Rd(iainia)/2.
1449
1450	! We assume that cs_star equals ci_star (see Comment in legend of Fig. 6 of Yin et al. (2009)
1451	! Equation 26 of Yin et al. (2009)
1452	Cs_star = (gbs(jv) * low_gamma_star(iainia) * Oi - ( 1. + &
1453	low_gamma_star(iainia) * alpha(jv) / 0.047) * Rd(iainia) + Rm ) &
1454	/ ( gbs(jv) + kp(jv) )
1455
1456
1457	! eq. 11 of Yin et al (2009)
1458	J2 = JJ(iainia) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) )
1459
1460	! Equation right after eq. (20d) of Yin et al. (2009)
1461	z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc ))
1462
1463	VpJ2 = fpsir(jv) * J2 * z / 2.
1464
1465	A_3=9999.
1466
1467	! See eq. right after eq. 18 of Yin et al. (2009)
1468	! The loop over limit_photo calculates the assimultion
1469	! rate of the Rubisco-limited CsO2 assimilation (Ac) and
1470	! the rate of electron transport-limited CO2 assimultion
1471	! (Aj). The overall assimilation will be the slower
1472	! of these two processes.
1473	DO limit_photo=1,2
1474	! Is Vc limiting the Assimilation
1475	IF ( limit_photo .EQ. 1 ) THEN
1476	a = 1. + kp(jv) / gbs(jv)
1477	b = 0.
1478	x1 = Vc2(iainia)
1479	x2 = KmC(iainia)/KmO(iainia)
1480	x3 = KmC(iainia)
1481	! Is J limiting the Assimilation
1482	ELSE
1483	a = 1.
1484	b = VpJ2
1485	x1 = (1.- fpsir(jv)) * J2 * z / 3.
1486	x2 = 7. * low_gamma_star(iainia) / 3.
1487	x3 = 0.
1488	ENDIF
1489
1490	m=fvpd(iainia)-g0(jv)/gb(iainia)
1491	d=g0(jv)(Ca(iainia)-Cs_star) + fvpd(iainia)Rd(iainia)
1492	f=(b-Rm-low_gamma_star(iainia)Oigbs(jv))x1d + agbs(jv)x1Ca(iainia)d
1493	j=(b-Rm+gbs(jv)x3 + x2gbs(jv)Oi)m + (alpha(jv)x2/0.047-1.)d &
1494	+ agbs(jv)(Ca(iainia)*m - d/gb(iainia) - (Ca(iainia) - Cs_star ))
1495
1496	g=(b-Rm-low_gamma_star(iainia)Oigbs(jv))x1m - (alpha(jv)low_gamma_star(iainia)/0.047+1.)x1*d &
1497	+ agbs(jv)x1(Ca(iainia)m - d/gb(iainia) - (Ca(iainia)-Cs_star ))
1498
1499	h=-((alpha(jv)low_gamma_star(iainia)/0.047+1.)x1m + (agbs(jv)x1(m-1.))/gb(iainia) )
1500	i= ( b-Rm + gbs(jv)x3 + x2gbs(jv)Oi )d + agbs(jv)Ca(iainia)*d
1501	l= ( alpha(jv)x2/0.047 - 1.)m - (agbs(jv)(m-1.))/gb(iainia)
1502
1503	p = (j-(h-l*Rd(iainia))) / l
1504	q = (i+j*Rd(iainia)-g) / l
1505	r = -(f-i*Rd(iainia)) / l
1506
1507	! See Yin et al. (2009) and Baldocchi (1994)
1508	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
1509	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
1510
1511	IF (QQ .GE. 0._r_std) THEN
1512	IF (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) THEN
1513	PSI = ACOS(UU/(QQ**1.5_r_std))
1514	A_3_tmp = -2._r_std * SQRT(QQ) * COS(( PSI + 4._r_std * PI)/3._r_std ) - p / 3._r_std
1515	IF (( A_3_tmp .LT. A_3 )) THEN
1516	A_3 = A_3_tmp
1517	ELSE
1518	! In case, J is not limiting the assimilation
1519	! we have to re-initialise a, b, x1, x2 and x3 values
1520	! in agreement with a Vc-limited assimilation
1521	a = 1. + kp(jv) / gbs(jv)
1522	b = 0.
1523	x1 = Vc2(iainia)
1524	x2 = KmC(iainia)/KmO(iainia)
1525	x3 = KmC(iainia)
1526	ENDIF
1527	ENDIF
1528	ELSE
1529	WRITE(numout,*) 'WARNING: unexpected value for QQ (sechiba_hydrol_arch.f90)'
1530	ENDIF
1531
1532	IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(iainia)) ) ) THEN
1533	WRITE(numout,*) 'We have a problem in jdiffuco_trans_co2'
1534	WRITE(numout,*) 'no real positive solution found for pft:',jv
1535	WRITE(numout,*) 'temp_sol:',temp_sol(iainia)
1536	WRITE(numout,*) 'vpd:',vpd(iainia)
1537	A_3 = -Rd(iainia)
1538	ENDIF
1539	assimi(iainia) = A_3
1540
1541	!++++ TEMP +++++
1542	IF(( x1 - (assimi(iainia) + Rd(iainia) )) == zero)THEN
1543
1544	! This causes a problem, since we really should not
1545	! be dividing by zero. However, it seems to happen
1546	! fairly rarely, so all we'll do is cycle to the
1547	! next loop and keep track of this warning so we
1548	! can look at it afterwards.
1549	warnings(iainia,jv,iwphoto)=warnings(iainia,jv,iwphoto)+un
1550
1551	IF(ld_warn)THEN
1552	WRITE(numout,*) 'jdiffuco, Getting ready to divide by zero! C4'
1553	WRITE(numout,*) '( x1 - ( assimi(iainia) + Rd(iainia) ) ', &
1554	( x1 - (assimi(iainia) + Rd(iainia) ))
1555	WRITE(numout,*) 'x1,assimi(iainia),Rd(iainia): ', &
1556	x1,assimi(iainia),Rd(iainia)
1557	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
1558	ENDIF
1559
1560	CYCLE
1561
1562	ENDIF
1563
1564	! Eq. 24 of Yin et al. (2009)
1565	Obs = ( alpha(jv) * assimi(iainia) ) / ( 0.047 * gbs(jv) ) + Oi
1566	! Eq. 23 of Yin et al. (2009)
1567	Cc = ( ( assimi(iainia) + Rd(iainia) ) * ( x2 * Obs + x3 ) + low_gamma_star(iainia) &
1568	* Obs * x1 ) &
1569	/ ( x1 - ( assimi(iainia) + Rd(iainia) ))
1570	! Eq. 22 of Yin et al. (2009)
1571	new_leaf_ci(iainia,jv,ilev) = ( Cc - ( b - assimi(iainia) - Rm ) / gbs(jv) ) / a
1572
1573
1574
1575	! JR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1576
1577	! For the cases in which the supply term condition applies, we need
1578	! to work backwards and calculate gs, and hence assimi for the grid
1579	! square and PFT in question, from the value of vbeta3 derived in the
1580	! enerbil code.
1581
1582	! This sequence is for the c4 vegetation types, and should be identical
1583	! to the c3 loop below, except for the definition of assimi
1584
1585	IF (hydrol_flag3(iainia, jv)) THEN
1586
1587	!WRITE(numout,*) '231113 in c4 loop ilev is: ', ilev, 'PFT is: ', jv
1588
1589	! WRITE(numout,*) '151113 in jenerbil_trans_co2 hydrol_flag3 C4 plants'
1590
1591	speed = MAX(min_wind, jwind(iainia))
1592
1593	! first case
1594
1595	IF (zqsvegrap(iainia) .lt. 1.0d0) THEN
1596
1597	!cresist = vbeta3(iainia, jv) / &
1598	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) + &
1599	! & veget_max(iainia, jv) * zqsvegrap(iainia) )
1600
1601	cresist = vbeta3(iainia, jv) / &
1602	& (veget_max(iainia, jv) * ( un - zqsvegrap(iainia)) + &
1603	& veget(iainia,jv) * zqsvegrap(iainia) )
1604
1605
1606	! WRITE(numout,*) 'cresist condition a)'
1607
1608	ELSE
1609
1610	! cresist = vbeta3(iainia, jv) / &
1611	! & ( veget_max(iainia, jv) * zqsvegrap(iainia) )
1612
1613	cresist = 0.0d0
1614
1615	! WRITE(numout,*) 'cresist condition b)'
1616
1617	END IF
1618
1619	IF (vbeta23(iainia, jv) .le. &
1620	& (veget(iainia, jv) * zqsvegrap(iainia) * cresist) ) THEN
1621	! do nothing
1622	ELSE
1623
1624	IF (zqsvegrap(iainia) .lt. 1.0d0) THEN
1625	!cresist = (vbeta3(iainia, jv) - vbeta23(iainia, jv) ) / &
1626	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) )
1627
1628	cresist = vbeta3(iainia, jv) / &
1629	& (veget_max(iainia, jv) * ( un - zqsvegrap(iainia)) + &
1630	& veget(iainia,jv) * zqsvegrap(iainia) )
1631
1632
1633	! WRITE(numout,*) 'cresist condition c)'
1634
1635	ELSE
1636
1637	! cresist = vbeta3(iainia, jv) / &
1638	! & ( veget_max(iainia, jv) * zqsvegrap(iainia) )
1639
1640	cresist = 0.0d0
1641
1642	END IF
1643
1644	END IF !
1645
1646	!
1647	IF (cresist .ne. 0.0d0) THEN
1648	sum_rveget_rstruct(iainia,jv) = &
1649	& (1.0d0/cresist - 1.0d0) * (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
1650	ELSE
1651	! sum_rveget_rstruct(iainia,jv) = &
1652	! & (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
1653
1654	sum_rveget_rstruct(iainia,jv) = 0.0d0
1655
1656	END IF
1657
1658	!WRITE(numout,*) '241113 sub_rveget_rstruct is: ', sum_rveget_rstruct(iainia,jv)
1659
1660	vbeta3pot(iainia,jv) = MAX(zero, veget_max(iainia,jv) * cresist)
1661
1662
1663	! calculate approximate gs, working backwards from vbeta3
1664	IF (sum_rveget_rstruct(iainia,jv) .ne. 0.0d0) THEN
1665	gstot(iainia) = 1.0d0 / sum_rveget_rstruct(iainia,jv)
1666	ELSE
1667	! Need to convert the units of g0
1668	gstot(iainia) = g0(jv)mmol_to_m_1 (temp_sol(iainia)/tp_00)*&
1669	(pb_std/pb(iainia))*ratio_H2O_to_CO2
1670	END IF
1671
1672	!WRITE(numout,*) '241113 gstot is: ', gstot
1673
1674	IF (gstot(iainia) .gt. 0.0d0) THEN
1675	! nothing
1676	ELSE
1677	!WRITE(,) '121113 gstot .le. 0'
1678	ENDIF
1679
1680
1681	! we now do a unit conversion, from m/s to mol/m^2/s, the reverse
1682	! of the calculation of section 2.5, below
1683
1684	gstot(iainia) = gstot(iainia) / (mmol_to_m_1 (temp_sol(iainia)/tp_00) &
1685	(pb_std/pb(iainia))*ratio_H2O_to_CO2)
1686
1687
1688	! gs(iainia) = gstot(iainia) / lai(iainia, jv)
1689	! gs(iainia) = gstot(iainia) * gs_distribution(iainia, jv, ilev) !/ lai_total(iainia,jv)
1690
1691	!**********************CHECK KIM: gs should not be lower than g0**********************************
1692
1693
1694	gs(iainia) = gstot(iainia)*gstot_frac(iainia,jv,ilev) / lai_per_level(iainia,jv,ilev)
1695
1696	! IF (gs(iainia) .gt. gstot(iainia)) THEN
1697	! gs(iainia) = gstot(iainia)
1698	! END IF
1699
1700
1701	IF (gs(iainia) .GE. g0(jv)) THEN
1702	assimi(iainia) = (((gs(iainia) - g0(jv)) * (Ca(iainia) - Cs_star)) / fvpd(iainia)) &
1703	& - Rd(iainia)
1704	ELSE
1705	assimi(iainia) = -Rd(iainia)
1706	ENDIF
1707
1708	assimi_store(iainia, jv, ilev) = assimi(iainia)
1709
1710	! WRITE(numout,*) '281113 assimi is: ', assimi
1711
1712	IF (assimi(iainia) .lt. min_sechiba) THEN
1713	assimi(iainia) = 0.0d0
1714	! WRITE(numout,) 'assimi(iainia) is negative **************************************'
1715	! WRITE(numout,*) 'assimi(iainia) is: ', assimi(iainia)
1716	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
1717	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
1718	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
1719	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
1720	! WRITE(numout,*) 'Cs_star is: ', Cs_star
1721	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
1722	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
1723	! WRITE(numout,) '********************************************************'
1724	ELSE IF (assimi(iainia) .gt. 500.0d0) THEN
1725	! WRITE(numout,*) 'assimi(iainia) is very large'
1726	! WRITE(numout,*) 'assimi(iainia) is: ', assimi(iainia)
1727	! WRITE(numout,*) 'temp_sol(iainia) is: ', temp_sol(iainia)
1728	! WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
1729	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
1730	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
1731	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
1732	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
1733	! WRITE(numout,*) 'Cs_star is: ', Cs_star
1734	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
1735	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
1736	! WRITE(numout,*) 'gstot(iainia) is: ', gstot(iainia)
1737	! WRITE(numout,*) 'sum_rveget_rstruct(iainia,jv) is: ', sum_rveget_rstruct(iainia,jv)
1738	! WRITE(numout,*) 'speed is: ', speed
1739	! WRITE(numout,*) 'q_cdrag(iainia) is: ',q_cdrag(iainia)
1740	! WRITE(numout,*) 'rveg_pft(jv) is: ', rveg_pft(jv)
1741	! WRITE(numout,*) 'gstot_frac(iainia,jv,ilev) is: ', gstot_frac(iainia,jv,ilev)
1742	! WRITE(numout,*) 'lai_per_level(iainia,jv,ilev) is: ', lai_per_level(iainia,jv,ilev)
1743	! WRITE(numout,) '********************************************************'
1744	END IF
1745
1746	ENDIF ! (hydrol_flag3(iainia, jv))
1747
1748	! with new assimilation rates for the grid squares and the PFTs for
1749	! which the supply term restriction replies, we let the rest of the
1750	! routine run its course from here, to calculate gs and gpp (vbeta3
1751	! remains the same as that calculated within enerbil)
1752
1753	! JR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1754
1755
1756	! Eq. 25 of Yin et al. (2009)
1757	! It should be Cs instead of Ca but it seems that
1758	! other equations in Appendix C make use of Ca
1759	! IF ( ABS( assimi(iainia) + Rd(iainia) ) .LT. min_sechiba ) THEN
1760
1761	!************CHECK KIM: shouldn't this be done before the calculation of assimi? *********************************
1762	IF ( gs(iainia) .LT. min_sechiba ) THEN
1763	gs(iainia) = g0(jv)
1764	ENDIF
1765	!
1766
1767	ENDDO
1768	!
1769	ELSE
1770	!
1771	! @addtogroup Photosynthesis
1772	! @{
1773	!
1774	! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
1775	!! \latexonly
1776	!! \input{jdiffuco_trans_co2_2.4.3.tex}
1777	!! \endlatexonly
1778	! @}
1779	!
1780	!
1781	!
1782	A_1=9999.
1783
1784	! See eq. right after eq. 18 of Yin et al. (2009)
1785	DO limit_photo=1,2
1786	! Is Vc limiting the Assimilation
1787	IF ( limit_photo .EQ. 1 ) THEN
1788	x1 = vc2(iainia)
1789	! It should be O not Oi (comment from Vuichard)
1790	x2 = KmC(iainia) * ( 1. + Oi / KmO(iainia) )
1791	! Is J limiting the Assimilation
1792	ELSE
1793	x1 = JJ(iainia)/4.
1794	x2 = 2. * gamma_star(iainia)
1795	ENDIF
1796
1797
1798	! See Appendix B of Yin et al. (2009)
1799	a = g0(jv) * ( x2 + gamma_star(iainia) ) + ( fvpd(iainia) ) * ( x1 - Rd(iainia) )
1800	b = Ca(iainia) * ( x1 - Rd(iainia) ) - gamma_star(iainia) * x1 - Rd(iainia) * x2
1801	c = Ca(iainia) + x2 + ( 1./gb(iainia) ) * ( x1 - Rd(iainia) )
1802	d = x2 + gamma_star(iainia)
1803	m = fvpd(iainia) / gb(iainia)
1804
1805	p = - ( d + a / gb(iainia) + fvpd(iainia) * c ) / m
1806
1807	q = ( d * ( x1 - Rd(iainia) ) + ac + ( fvpd(iainia) ) b ) / m
1808	r = - a * b / m
1809
1810	! See Yin et al. (2009)
1811	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
1812	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
1813
1814	IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN
1815	PSI = ACOS(UU/(QQ**1.5_r_std))
1816	A_1_tmp = -2._r_std * SQRT(QQ) * COS( PSI / 3._r_std ) - p / 3._r_std
1817	IF (( A_1_tmp .LT. A_1 )) THEN
1818	A_1 = A_1_tmp
1819	ELSE
1820	! In case, J is not limiting the assimilation
1821	! we have to re-initialise x1 and x2 values
1822	! in agreement with a Vc-limited assimilation
1823	x1 = vc2(iainia)
1824	! It should be O not Oi (comment from Vuichard)
1825	x2 = KmC(iainia) * ( 1. + Oi / KmO(iainia) )
1826	ENDIF
1827	ENDIF
1828	ENDDO
1829
1830	IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(iainia)) ) ) THEN
1831	WRITE(numout,*) 'We have a problem in jdiffuco_trans_co2'
1832	WRITE(numout,*) 'no real positive solution found for pft:',jv
1833	WRITE(numout,*) 'temp_sol:',temp_sol(iainia)
1834	WRITE(numout,*) 'vpd:',vpd(iainia)
1835	WRITE(numout,*) 'Setting the solution to -Rd(iainia):',-Rd(iainia)
1836	A_1 = -Rd(iainia)
1837	ENDIF
1838	assimi(iainia) = A_1
1839	! Eq. 18 of Yin et al. (2009)
1840
1841	!++++++++ TEMP +++++++
1842	IF(( x1 - (assimi(iainia) + Rd(iainia) )) == zero)THEN
1843	! This causes a problem, since we really should not
1844	! be dividing by zero. However, it seems to happen
1845	! fairly rarely, so all we'll do is cycle to the
1846	! next loop and keep track of this warning so we
1847	! can look at it afterwards.
1848	warnings(iainia,jv,iwphoto)=warnings(iainia,jv,iwphoto)+un
1849
1850	IF(ld_warn)THEN
1851	WRITE(numout,*) 'jdiffuco, Getting ready to divide by zero! C3'
1852	WRITE(numout,*) '( x1 - ( assimi(iainia) + Rd(iainia) ) ', &
1853	( x1 - (assimi(iainia) + Rd(iainia) ))
1854	WRITE(numout,*) 'x1,assimi(iainia),Rd(iainia): ', &
1855	x1,assimi(iainia),Rd(iainia)
1856	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
1857	ENDIF
1858
1859	CYCLE
1860
1861	ENDIF
1862	!++++++++++++++++++++++++++++++
1863
1864	Cc = ( gamma_star(iainia) * x1 + ( assimi(iainia) + Rd(iainia) ) * x2 ) &
1865	/ ( x1 - ( assimi(iainia) + Rd(iainia) ) )
1866	! Eq. 17 of Yin et al. (2009)
1867	! new_leaf_ci(iainia,jv,ilev) = Cc
1868	new_leaf_ci(iainia,jv,ilev) = Cc
1869	! See eq. right after eq. 15 of Yin et al. (2009)
1870	ci_star = gamma_star(iainia)
1871	!
1872
1873
1874	! JR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1875
1876	! For the cases in which the supply term condition applies, we need
1877	! to work backwards and calculate gs, and hence assimi for the grid
1878	! square and PFT in question, from the value of vbeta3 derived in the
1879	! enerbil code.
1880
1881	! This sequence is for the c3 vegetation types, and should be identical
1882	! to the c4 loop above, except for the definition of assimi
1883
1884	IF (hydrol_flag3(iainia, jv)) THEN
1885
1886	!WRITE(numout,*) 'in jdiffuco_trans_co2'
1887	!WRITE(numout,*) '231113 in c3 loop ilev is: ', ilev, 'PFT is: ', jv
1888
1889	speed = MAX(min_wind, jwind(iainia))
1890
1891	! first case
1892	! WRITE(numout,*) '151113 cresist before: ', cresist
1893
1894	IF (zqsvegrap(iainia) .lt. 1.0d0) THEN
1895	!calculation of cresist may need to be check
1896	! Here may to to devided gradient of spacific
1897	! humidity, at least which need to consist with the
1898	! caculation in the diffuco_trans_co2..
1899	!cresist = vbeta3(iainia, jv) / &
1900	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) + &
1901	! & veget_max(iainia, jv) * zqsvegrap(iainia) )
1902	cresist = vbeta3(iainia, jv) / &
1903	& (veget_max(iainia, jv) * ( un - zqsvegrap(iainia)) + &
1904	& veget(iainia,jv) * zqsvegrap(iainia) )
1905
1906
1907	! WRITE(numout,*) 'cresist condition a)'
1908
1909	ELSE
1910
1911	! cresist = vbeta3(iainia, jv) / &
1912	! & ( veget_max(iainia, jv) * zqsvegrap(iainia) )
1913
1914	cresist = 0.0d0
1915
1916	! WRITE(numout,*) 'cresist condition b)'
1917
1918	END IF
1919
1920
1921	IF (vbeta23(iainia, jv) .le. &
1922	& (veget(iainia, jv) * zqsvegrap(iainia) * cresist) ) THEN
1923
1924	! do nothing
1925
1926	ELSE
1927
1928	IF (zqsvegrap(iainia) .lt. 1.0d0) THEN
1929
1930	!cresist = (vbeta3(iainia, jv) - vbeta23(iainia, jv) ) / &
1931	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) )
1932
1933	cresist = vbeta3(iainia, jv) / &
1934	& (veget_max(iainia, jv) * ( un - zqsvegrap(iainia)) + &
1935	& veget(iainia,jv) * zqsvegrap(iainia) )
1936	! WRITE(numout,*) 'cresist condition c)'
1937
1938	ELSE
1939
1940	! cresist = vbeta3(iainia, jv) / &
1941	! & ( veget_max(iainia, jv) * zqsvegrap(iainia) )
1942
1943	cresist = 0.0d0
1944
1945	! WRITE(numout,*) 'cresist condition d)'
1946
1947	END IF
1948
1949	END IF !
1950
1951	!WRITE(numout,*) '241113 cresist is: ', cresist
1952
1953	!***************************TEMP************************
1954	!WRITE(,) '151113 cresist (jdiffuco_trans_co2) is: ', cresist
1955	!WRITE(,) '151113 speed: ', speed
1956	!WRITE(,) '151113 cresist (jdiffuco_trans_co2) is: ', q_cdrag(iainia)
1957	!WRITE(,) '151113 rveg_pft (jdiffuco_trans_co2) is: ', rveg_pft(jv)
1958	!***************************TEMP**************************
1959
1960	IF (cresist .ne. 0.0d0) THEN
1961	sum_rveget_rstruct(iainia,jv) = &
1962	& (1.0d0/cresist - 1.0d0) * (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
1963	ELSE
1964	sum_rveget_rstruct(iainia,jv) = &
1965	& (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
1966	END IF
1967
1968	IF (cresist .ne. 0.0d0) THEN
1969	sum_rveget_rstruct(iainia,jv) = &
1970	& (1.0d0/cresist - 1.0d0) * (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
1971	ELSE
1972	! sum_rveget_rstruct(iainia,jv) = &
1973	! & (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
1974
1975	sum_rveget_rstruct(iainia,jv) = 0.0d0
1976
1977	END IF
1978
1979	!WRITE(numout,*) '241113 sub_rveget_rstruct is: ', sum_rveget_rstruct(iainia,jv)
1980
1981	vbeta3pot(iainia,jv) = MAX(zero, veget_max(iainia,jv) * cresist)
1982
1983
1984	! calculate approximate gs, working backwards from vbeta3
1985	!
1986	IF (sum_rveget_rstruct(iainia,jv) .ne. 0.0d0) THEN
1987	gstot(iainia) = 1.0d0 / sum_rveget_rstruct(iainia,jv)
1988	ELSE
1989	! Need to convert the units of g0
1990	gstot(iainia) = g0(jv)mmol_to_m_1 (temp_sol(iainia)/tp_00)*&
1991	(pb_std/pb(iainia))*ratio_H2O_to_CO2
1992	END IF
1993
1994
1995
1996	!******CHECK KIM: this what I had for a divide by zero bug fix, for now I keep what's in James update
1997	! IF (cresist .NE. un .AND. speed .GT. min_sechiba .AND. &
1998	! & q_cdrag(iainia) .GT. min_sechiba .AND. rveg_pft(jv) .GT. min_sechiba) THEN
1999	!
2000	! sum_rveget_rstruct(iainia,jv) = &
2001	! & (1.0d0/cresist - 1.0d0) * (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
2002	!ENDIF
2003
2004	!WRITE(numout,*) '241113 sub_rveget_rstruct is: ', sum_rveget_rstruct(iainia,jv)
2005
2006
2007	IF (cresist .lt. min_sechiba) THEN
2008	IF(ld_warn .OR. ld_wstress) THEN
2009	! WRITE(numout,) 'cresist is very small or negative *********************'
2010	! WRITE(numout,*) 'cresist is: ', cresist
2011	! WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
2012	! WRITE(numout,*) 'veget_max(iainia, jv) is: ', veget_max(iainia, jv)
2013	! WRITE(numout,*) 'zqsvegrap(iainia) is: ', zqsvegrap(iainia)
2014	! WRITE(numout,) '********************************************************'
2015	ENDIF
2016	ELSE IF (assimi(iainia) .gt. 100.0d0) THEN
2017	WRITE(numout,*) 'assimi(iainia) is very large'
2018	END IF
2019
2020
2021	IF (gstot(iainia) .gt. 0.0d0) THEN
2022	! nothing
2023	ELSE
2024	!WRITE(,) '121113 gstot .le. 0'
2025	ENDIF
2026
2027	! we assume lai is the sum of lai_per_level over all of the levels
2028
2029	! we now do a unit conversion, from m/s to mol/m^2/s, the reverse
2030	! of the calculation of section 2.5, below
2031
2032	gstot(iainia) = gstot(iainia) / (mmol_to_m_1 (temp_sol(iainia)/tp_00) &
2033	(pb_std/pb(iainia))*ratio_H2O_to_CO2)
2034
2035
2036
2037
2038	gs(iainia) = gstot(iainia)*gstot_frac(iainia,jv,ilev) / lai_per_level(iainia,jv,ilev)
2039
2040	IF (gs(iainia) .gt. gstot(iainia)) THEN
2041	gs(iainia) = gstot(iainia)
2042	END IF
2043
2044
2045
2046	! define gs, based on vbeta3
2047
2048	IF (gs(iainia) .GE. g0(jv)) THEN
2049	assimi(iainia) = (((gs(iainia) - g0(jv)) * (new_leaf_ci(iainia,jv,ilev) - ci_star)) / fvpd(iainia)) &
2050	& - Rd(iainia)
2051	ELSE
2052	assimi(iainia) = -Rd(iainia)
2053
2054	ENDIF
2055
2056	! assimi(iainia) = (gs(iainia) - g0(jv)) * (new_leaf_ci(iainia,jv,ilev) - ci_star) / fvpd(iainia) &
2057	! & - Rd(iainia)
2058
2059	assimi_store(iainia, jv, ilev) = assimi(iainia)
2060
2061
2062	! IF (assimi(iainia) .lt. min_sechiba) THEN
2063	! WRITE(numout,*) '281113 assimi is: ', assimi
2064	! END IF
2065
2066	IF (assimi(iainia) .lt. min_sechiba) THEN
2067	assimi(iainia) = 0.0d0
2068	! WRITE(numout,) 'assimi(iainia) is negative **************************************'
2069	! WRITE(numout,*) 'assimi(iainia) is: ', assimi(iainia)
2070	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
2071	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
2072	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
2073	! WRITE(numout,*) 'ci_star is: ', ci_star
2074	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
2075	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
2076	! WRITE(numout,) '********************************************************'
2077	ELSE IF (assimi(iainia) .gt. 500.0d0) THEN
2078	! WRITE(numout,*) 'assimi(iainia) is very large'
2079	! WRITE(numout,*) 'assimi(iainia) is: ', assimi(iainia)
2080	! WRITE(numout,*) 'temp_sol(iainia) is: ', temp_sol(iainia)
2081	! WRITE(numout,*) 'vbeta3(iainia, jv) is: ', vbeta3(iainia, jv)
2082	! WRITE(numout,*) 'gs(iainia) is: ', gs(iainia)
2083	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
2084	! WRITE(numout,*) 'new_leaf_ci(iainia,jv,ilev) is: ', new_leaf_ci(iainia,jv,ilev)
2085	! WRITE(numout,*) 'ci_star is: ', ci_star
2086	! WRITE(numout,*) 'Ca(iainia) is: ', Ca(iainia)
2087	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
2088	! WRITE(numout,*) 'Rd(iainia) is: ', Rd(iainia)
2089	! WRITE(numout,*) 'gstot(iainia) is: ', gstot(iainia)
2090	! WRITE(numout,*) 'sum_rveget_rstruct(iainia,jv) is: ', sum_rveget_rstruct(iainia,jv)
2091	! WRITE(numout,*) 'speed is: ', speed
2092	! WRITE(numout,*) 'q_cdrag(iainia) is: ',q_cdrag(iainia)
2093	! WRITE(numout,*) 'rveg_pft(jv) is: ', rveg_pft(jv)
2094	! WRITE(numout,*) 'gstot_frac(iainia,jv,ilev) is: ', gstot_frac(iainia,jv,ilev)
2095	! WRITE(numout,*) 'lai_per_level(iainia,jv,ilev) is: ', lai_per_level(iainia,jv,ilev)
2096	! WRITE(numout,) '********************************************************'
2097	END IF
2098
2099
2100
2101
2102	! WRITE(numout,*) '151113 assimi(iainia) (jdiffuco_trans_co2) is: ', assimi(iainia)
2103
2104	ENDIF ! (hydrol_flag3(iainia, jv))
2105
2106	! with new assimilation rates for the grid squares and the PFTs for
2107	! which the supply term restriction replies, we let the rest of the
2108	! routine run its course from here, to calculate gs and gpp (vbeta3
2109	! remains the same as that calculated within enerbil)
2110
2111	! JR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2112
2113	! Eq. 25 of Yin et al. (2009)
2114	! It should be Cs instead of Ca but it seems that
2115	! other equations in Appendix C make use of Ca
2116	! IF ( ABS( assimi(iainia) + Rd(iainia) ) .LT. min_sechiba ) THEN
2117	IF ( gs(iainia) .LT. min_sechiba ) THEN
2118	gs(iainia) = g0(jv)
2119	ELSE
2120	! gs(iainia) = g0(jv) + ( assimi(iainia) + Rd(iainia) ) / ( Ca(iainia) - Cs_star ) * fvpd(iainia)
2121	gs(iainia) = gs(iainia) ! as calculated in loop
2122	ENDIF
2123
2124
2125
2126	ENDIF
2127
2128
2129	!
2130	! @addtogroup Photosynthesis
2131	! @{
2132	!
2133	!! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n
2134	!! \latexonly
2135	!! \input{jdiffuco_trans_co2_2.4.4.tex}
2136	!! \endlatexonly
2137	! @}
2138	!
2139	!
2140	! keep stomatal conductance of the levels which we decide
2141	! are the top levels. Notice that we multiply by the LAI
2142	! in the level here, so this is no longer a per leaf quantity,
2143	! but the whole quantity.
2144	IF ( top_level(iainia,ilev) ) THEN
2145	gs_top(iainia) = gs_top(iainia)+gs(iainia)*lai_per_level(iainia,jv,ilev)
2146	!
2147	ENDIF
2148	!
2149	! @addtogroup Photosynthesis
2150	! @{
2151	!
2152	!! 2.4.5 Integration at the canopy level\n
2153	!! \latexonly
2154	!! \input{jdiffuco_trans_co2_2.4.5.tex}
2155	!! \endlatexonly
2156	! @}
2157	! total assimilation and conductance
2158	!assimtot(iainia) = assimtot(iainia) + &
2159	! assimi(iainia) * (lai_per_level(ilev+1)-lai_per_level(ilev))
2160	!Rdtot(iainia) = Rdtot(iainia) + &
2161	! Rd(iainia) * (lai_per_level(ilev+1)-lai_per_level(ilev))
2162	!gstot(iainia) = gstot(iainia) + &
2163	! gs(iainia) * (lai_per_level(ilev+1)-lai_per_level(ilev))
2164	!
2165	!++++++++++++++++++++
2166	! The equations give the amount of carbon produced per m**2 if
2167	! leaf surface. Therefore, since we want GPP in units of m**2 of land surface,
2168	! we have to multiply by the LAI.
2169	assimtot(iainia) = assimtot(iainia) + &
2170	assimi(iainia) * lai_per_level(iainia,jv,ilev)
2171	! assimi(iainia)
2172	IF(test_pft == jv .AND. ld_photo)THEN
2173	WRITE(numout,*) 'assimi ',ilev,iainia,assimi(iainia),lai_per_level(iainia,jv,ilev)
2174	ENDIF
2175	Rdtot(iainia) = Rdtot(iainia) + &
2176	! Rd(iainia)
2177	Rd(iainia) * lai_per_level(iainia,jv,ilev)
2178	gstot(iainia) = gstot(iainia) + &
2179	gs(iainia) * lai_per_level(iainia,jv,ilev)
2180	! gs(iainia)
2181	! WRITE(numout,*) '151113 iainia is: ', iainia, 'gs(iainia) is: ', gs(iainia)
2182
2183	!
2184	gstot_output(iainia,jv,ilev) = gstot(iainia)
2185	gs_output(iainia,jv,ilev) = gs(iainia)
2186	!
2187
2188	ENDIF ! if there is absorbed light
2189	ENDDO ! loop over grid points
2190	ENDDO gridloop1 ! loop over LAI steps
2191
2192
2193
2194	IF (ld_warn .AND. assimtot(iainia) .gt. 20) THEN
2195	WRITE(numout,*) 'assimtot is very high'
2196	WRITE(numout,*) 'assimtot is: ', assimtot(iainia)
2197	DO ilev=1,nlevels_tot
2198	WRITE (numout,*) 'level: ', ilev, ' pft: ', jv , ' assimi_store: ', assimi_store(iainia, jv, ilev), &
2199	&' lai_per_level: ', lai_per_level(iainia, jv, ilev)
2200	END DO
2201	END IF
2202
2203
2204
2205
2206	DO jij = 1, kjpindex
2207	IF (hydrol_flag2(jij)) THEN
2208	! WRITE(,) '091113 new assimi(grid point) is: ', assimi(jij)
2209	END IF
2210	END DO !jij =1, nvm
2211
2212	! WRITE(numout,*) '151113b gstot is: ', gstot
2213	! WRITE(numout,*) '151113b lai_per_level is: ', lai_per_level
2214
2215	! Due to the water limited by the transpiration, we nned to recalculated
2216	! the stomata conductance(resistance) for the GPP.
2217	! In this section, the assimlation should be recalculated.
2218
2219	!! 2.5 Calculate resistances
2220
2221	gridloop2: DO inia=1,nia
2222	!
2223	iainia=index_assi(inia)
2224	!
2225	! We hit a fringe case above where we could not
2226	! calculate the GPP. Skip this loop so we don't
2227	! divide by zero.
2228	IF(lskip(iainia)) CYCLE
2229
2230	! DO inina=1,nina
2231	IF (hydrol_flag3(iainia, jv)) THEN
2232	!nothing
2233	ELSE
2234	CYCLE gridloop2
2235	END IF !(hydrol_flag3)
2236	! END DO !inina=1,nina
2237
2238
2239	!++++++ CHECK ++++++
2240	! We introduce ::coupled_lai to calculate the :: vbeta3
2241	! In an closed canopy not all the leaves fully interact with the
2242	! atmosphere because leaves can shelter each other. In a more open
2243	! canopy most leaves can interact with the atmosphere. We use the
2244	! ratio of lai over ::Pgap as a proxy for the amount of canopy that
2245	! can interact with the atmosphere. Clearly that amount of lai that
2246	! interacts should never exceed the:: total_lai.
2247	coupled_lai(iainia, jv) = biomass_to_coupled_lai( biomass(iainia, &
2248	jv, ileaf, icarbon), veget(iainia,jv), jv)
2249
2250	total_lai(iainia, jv) = biomass_to_lai( biomass(iainia, jv, ileaf,&
2251	icarbon), jv )
2252
2253	gstot(iainia) = gstot(iainia)* (coupled_lai(iainia,jv)/total_lai(iainia,jv))
2254
2255
2256	! WRITE(numout,*) '151113 inia is: ', inia, 'gstot is: ', gstot(iainia)
2257	!
2258	!! Mean stomatal conductance for CO2 (mol m-2 s-1)
2259	gsmean(iainia,jv) = gstot(iainia)
2260	!
2261	! cimean is the "mean ci" calculated in such a way that assimilation
2262	! calculated in enerbil is equivalent to assimtot
2263	!
2264	IF((assimtot(iainia)+Rdtot(iainia)) .GT. min_sechiba)THEN
2265	cimean(iainia,jv) = (gsmean(iainia,jv)-g0(jv)) &
2266	/ (fvpd(iainia)*(assimtot(iainia)+Rdtot(iainia))) &
2267	+ gamma_star(iainia)
2268	ELSE
2269	! According to Nicolas, this variable is completely diagnostic and not
2270	! used for anything. If the analytical photosynthesis is unable to find
2271	! a solution for the A_1 coefficients in any of the levels, assimtot will
2272	! be equal to -Rd, which means we will be dividing by zero. This loop
2273	! is to prevent that, and since cimean is diagnostic is doesn't matter
2274	! what the value is.
2275	cimean(iainia,jv) = val_exp
2276	ENDIF
2277
2278	! conversion from umol m-2 (PFT) s-1 to gC m-2 (mesh area) tstep-1
2279	gpp(iainia,jv) = assimtot(iainia)12e-6veget_max(iainia,jv)*dtradia
2280
2281	! IF (gpp(iainia, jv) .lt. min_sechiba) THEN
2282	! WRITE(numout,) 'gpp is very small **************************************'
2283	! WRITE(numout,*) 'assimtot(iainia) is: ', assimtot(iainia)
2284	! WRITE(numout,*) 'fvpd(iainia) is: ', fvpd(iainia)
2285	! WRITE(numout,*) 'Rdtot(iainia) is: ', Rdtot(iainia)
2286	! WRITE(numout,*) 'veget_max(iainia,jv) is: ', veget_max(iainia,jv)
2287	! WRITE(numout,*) 'gsmean(iainia,jv) is: ', gsmean(iainia,jv)
2288	! WRITE(numout,*) 'g0(jv) is: ', g0(jv)
2289	! WRITE(numout,) '********************************************************'
2290	! ELSE IF (gpp(iainia, jv) .gt. 100.0d0) THEN
2291	! WRITE(numout,*) 'gpp is very large'
2292	! END IF
2293
2294	!IF(test_grid ==1 .AND. jv == test_pft)THEN
2295	! istep=istep+1
2296	! WRITE(numout,*) 'gpp jdiffuco: ',istep, gpp(iainia,jv), assimtot(iainia),veget_max(iainia,jv)
2297	!ENDIF
2298
2299	! conversion from (mol m-2 s-1) to (umol m-2 s-1)
2300	gsmean(iainia,jv) = gsmean(iainia,jv)*1e-6
2301	!
2302	! conversion from mol/m^2/s to m/s
2303	!
2304	! As in Pearcy, Schulze and Zimmermann
2305	! Measurement of transpiration and leaf conductance
2306	! Chapter 8 of Plant Physiological Ecology
2307	! Field methods and instrumentation, 1991
2308	! Editors:
2309	!
2310	! Robert W. Pearcy,
2311	! James R. Ehleringer,
2312	! Harold A. Mooney,
2313	! Philip W. Rundel
2314	!
2315	! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online)
2316	!
2317	!
2318	!we have already do the unit convertoion in the assimilation loop (gridloop1)
2319	!
2320	!gstot(iainia) = mmol_to_m_1 (temp_sol(iainia)/tp_00)&
2321	! (pb_std/pb(iainia))gstot(iainia)ratio_H2O_to_CO2
2322	!gstop(iainia) = mmol_to_m_1 * (temp_sol(iainia)/tp_00)*&
2323	! (pb_std/pb(iainia))gs_top(iainia)ratio_H2O_to_CO2
2324
2325
2326	!+++++++++Output some information regarding ::rveget and pft
2327	!IF (jv == test_pft .AND. test_grid == 1) THEN
2328	! WRITE(numout, '(A20, ES15.5)' ) 'New_Vbeta3:', vbeta3(iainia,jv)
2329	! WRITE(numout, '(A20, I5)' ) 'pft type is: ', jv
2330	! WRITE(numout, '(A20, ES15.5)') 'New_Gs_total:', gstot(iainia)
2331	!ENDIF
2332	!+++++++++++++++++++++++++++++++++++
2333
2334	!rveget(iainia,jv) = un/gstop(iainia)
2335
2336	!+++++CHECK+++++++++++++++++++++++++
2337	!For the grass land, we assume that there are still
2338	!some stomata on the other side of its leaves.
2339	!IF (.NOT. is_tree(jv) ) THEN
2340	! gstot(iainia) = 1.25_r_std * gstot(iainia)
2341	!ENDIF
2342	!we try to set rvegt as un/gstot
2343	rveget(iainia,jv) = un/gstot(iainia)
2344
2345	IF (lai(iainia, jv) .lt. 0.1d0) THEN
2346	rveget_min(iainia,jv) = (defc_plus / kzero(jv)) * (un / 0.1d0)
2347	ELSE
2348	rveget_min(iainia,jv) = (defc_plus / kzero(jv)) * (un / lai(iainia,jv))
2349	END IF
2350
2351	!
2352	! rstruct is the difference between rtot (=1./gstot) and rveget
2353	!
2354	! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
2355	!rstruct(iainia,jv) = un/gstot(iainia) - &
2356	! rveget(iainia,jv)
2357	!rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
2358	! rveget(iainia,jv), min_sechiba)
2359	!
2360	!
2361	!! wind is a global variable of the diffuco module.
2362	speed = MAX(min_wind, jwind(iainia))
2363	!
2364	! beta for transpiration
2365	!
2366	! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
2367	!! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
2368	!! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
2369	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2370	! (un - zqsvegrap(iainia)) * &
2371	! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
2372	! rstruct(iainia,jv))))
2373
2374
2375	!! Global resistance of the canopy to evaporation
2376	cresist=(un / (un + speed * q_cdrag(iainia) * &
2377	(rveg_pft(jv)*(rveget(iainia,jv) ))))
2378
2379	!WRITE(,) '051113 within jdiffuco_trans_co2 before defn. vbeta3(1,:) is', vbeta3(1,:)
2380
2381	! 051113 commenting out the calculation of vbeta3 (as this is already defined, either from
2382	! the first call to this module or as part of the new enerbil (when constrained by the supply
2383	! term))
2384
2385	! vbeta3(icinic,jv) = veget_max(icinic,jv) * &
2386	! (un - zqsvegrap(icinic)) * cresist + &
2387	! !! $ ! Addition Nathalie - June 2006
2388	! !! $ vbeta3(icinic,jv) = vbeta3(icinic,jv) + &
2389	! ! Corrections Nathalie - 09 November 2009 : veget => veget_max
2390	! ! MIN( vbeta23(icinic,jv), veget(icinic,jv) * &
2391	! MIN( vbeta23(icinic,jv), veget_max(icinic,jv) * &
2392	! ! zqsvegrap(icinic) * humrel(icinic,jv) * &
2393	! zqsvegrap(icinic) * cresist )
2394	! Fin ajout Nathalie
2395
2396	!WRITE(,) '051113 within jdiffuco_trans_co2 after defn. vbeta3(1,:) is', vbeta3(1,:)
2397
2398
2399
2400
2401	! vbeta3pot for computation of potential transpiration (needed for irrigation)
2402	! vbeta3pot(iainia,jv) = MAX(zero, veget_max(iainia,jv) * cresist)
2403
2404
2405
2406
2407	! 211013 - we have derived a new vbeta3 from the enerbil segment
2408
2409
2410	! WRITE(,) '051113 2nd call: vbeta3(iainia, jv) is: ', vbeta3(icinic, jv)
2411	! WRITE(,) '051113 2nd call: veget_max(iainia, jv) is: ', veget_max(icinic, jv)
2412	! WRITE(,) '051113 2nd call: zqsvegrap(iainia) is: ', zqsvegrap(icinic)
2413	! WRITE(,) '051113 2nd call: vbeta23(iainia, jv) is: ', vbeta23(icinic, jv)
2414	! WRITE(,) '051113 2nd call: cresist is: ', cresist
2415
2416
2417
2418
2419
2420
2421	! first case
2422	!WRITE(,) '051113 previous cresist: ', cresist
2423	!cresist = vbeta3(iainia, jv) / &
2424	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) + &
2425	! & veget_max(iainia, jv) * zqsvegrap(iainia) )
2426
2427
2428	!IF (vbeta23(iainia, jv) .le. &
2429	! & (veget_max(iainia, jv) * zqsvegrap(iainia) * cresist) ) THEN
2430	! do nothing
2431	! ELSE
2432	! cresist = vbeta3(iainia, jv) / &
2433	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) + &
2434	! & vbeta23(iainia, jv) )
2435	!cresist = (vbeta3(iainia, jv) - vbeta23(iainia, jv) ) / &
2436	! & ( veget_max(iainia, jv) * (1.0d0 - zqsvegrap(iainia)) )
2437
2438	!END IF !
2439
2440
2441
2442
2443
2444
2445
2446	!WRITE(,) '051113 new cresist: ', cresist
2447
2448
2449	!! Potentiometer to set vegetation resistance (unitless)
2450
2451	!sum_rveget_rstruct(iainia,jv) = &
2452	! & (1.0d0/cresist) * (1.0d0 / (speed * q_cdrag(iainia) * rveg_pft(jv)))
2453
2454
2455
2456
2457
2458
2459
2460
2461	!WRITE(,) '051113 previous vbeta3pot: ', vbeta3pot
2462	!vbeta3pot(iainia,jv) = MAX(zero, veget_max(iainia,jv) * cresist)
2463	!WRITE(,) '051113 new vbeta3pot: ', vbeta3pot
2464
2465
2466
2467
2468
2469
2470
2471	! vbetaco2 is not used anymore, replaced by gsmean, so I need to rethink the calculation
2472	! of these co-efficients when restriction by the supply term applies
2473
2474	! WRITE(,) '211013c previous vbetaco2: ', vbetaco2
2475	! vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
2476	! (un / (un + speed * q_cdrag(iainia) * &
2477	! (sum_rveget_rstruct(iainia,jv)))) / O2toCO2_stoechio
2478	! WRITE(,) '211013c new vbetaco2: ', vbetaco2
2479
2480	! cimean(iainia,jv) = ccanopy(iainia) - &
2481	! assimtot(iainia) / &
2482	! ( vbetaco2(iainia,jv)/veget_max(iainia,jv) * &
2483	! rau(iainia) * speed * q_cdrag(iainia))
2484
2485
2486
2487
2488	! ! beta for transpiration
2489	! !
2490	! ! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
2491	! !! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
2492	! !! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
2493	! !vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2494	! ! (un - zqsvegrap(iainia)) * &
2495	! ! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
2496	! ! rstruct(iainia,jv))))
2497	! !! Global resistance of the canopy to evaporation
2498	! cresist=(un / (un + speed * q_cdrag(iainia) * &
2499	! (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
2500
2501
2502	! vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2503	! (un - zqsvegrap(iainia)) * cresist + &
2504	!!!$ ! Addition Nathalie - June 2006
2505	!!!$ vbeta3(iainia,jv) = vbeta3(iainia,jv) + &
2506	! ! Corrections Nathalie - 09 November 2009 : veget => veget_max
2507	!! MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
2508	! MIN( vbeta23(iainia,jv), veget_max(iainia,jv) * &
2509	!! zqsvegrap(iainia) * humrel(iainia,jv) * &
2510	! zqsvegrap(iainia) * cresist )
2511	! ! Fin ajout Nathalie
2512
2513	! ! vbeta3pot for computation of potential transpiration (needed for irrigation)
2514	! vbeta3pot(iainia,jv) = MAX(zero, veget_max(iainia,jv) * cresist)
2515	!
2516	! beta for assimilation. The difference is that surface
2517	! covered by rain (un - zqsvegrap(iainia)) is not taken into account
2518	! The factor 1.6 is conversion for H2O to CO2 conductance.
2519	! vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
2520	! (un / (un + q_cdrag(iainia) * &
2521	! (rveget(iainia,jv))))/1.6
2522	!
2523
2524	! !! Global resistance of the canopy to CO2 transfer
2525	! vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
2526	! (un / (un + speed * q_cdrag(iainia) * &
2527	! (rveget(iainia,jv) + rstruct(iainia,jv)))) / O2toCO2_stoechio
2528	!
2529	! cimean is the "mean ci" calculated in such a way that assimilation
2530	! calculated in enerbil is equivalent to assimtot
2531	!
2532	! cimean(iainia,jv) = ccanopy(iainia) - &
2533	! assimtot(iainia) / &
2534	! ( vbetaco2(iainia,jv)/veget_max(iainia,jv) * &
2535	! rau(iainia) * speed * q_cdrag(iainia))
2536
2537	!
2538	!
2539	ENDDO gridloop2 !loop over grid points
2540	!
2541	ENDIF ! if nia LT zero
2542	!
2543	ENDDO ! loop over vegetation types
2544	!
2545	IF (long_print) WRITE (numout,*) ' jdiffuco_trans_co2 done '
2546
2547	!WRITE (numout,*) 'calling jnetcdf_trans_co2'
2548
2549	! this writes canopy profile data from this routine to a seperate netCDF file every time step,
2550	! though it does slow the runtime considerably, so is only for diagnostics
2551
2552	!CALL jnetcdf_trans_co2(kjpindex, james_netcdf_flag, gs_distribution(1,:,:), gs_output(1,:,:), &
2553	! & lai_per_level(1,:,:), light_tran_to_level(1,:,:), gs_diffuco_output(1,:,:), &
2554	! & gstot_output(1,:,:))
2555
2556
2557
2558
2559	IF(ld_photo)THEN
2560	WRITE(numout,*) 'gpp at the end of jdiffuco_trans_co2 ',gpp
2561	ENDIF
2562
2563
2564	END SUBROUTINE jdiffuco_trans_co2
2565
2566
2567
2568	! ---------------------------------------------------------------------------------------
2569
2570
2571
2572	FUNCTION Arrhenius (kjpindex,temp,ref_temp,energy_act) RESULT ( val_arrhenius )
2573	!! 0.1 Input variables
2574
2575	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
2576	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
2577	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
2578	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
2579
2580	!! 0.2 Result
2581
2582	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
2583	!! a Arrhenius function (-)
2584
2585	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:))))
2586	END FUNCTION Arrhenius
2587
2588	FUNCTION Arrhenius_modified (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
2589	!! 0.1 Input variables
2590
2591	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
2592	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
2593	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
2594	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
2595	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
2596	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
2597
2598	!! 0.2 Result
2599
2600	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
2601	!! a Arrhenius function (-)
2602
2603	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
2604	* (1. + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * RR ))) &
2605	/ (1. + EXP( (temp(:) * entropy(:) - energy_deact) / ( RR*temp(:))))
2606
2607	END FUNCTION Arrhenius_modified
2608
2609
2610	! ---------------------------------------------------------------------------------------
2611
2612
2613	SUBROUTINE jnetcdf_trans_co2(kjpindex, james_netcdf_flag, gs_distribution_in, gs_output_in, &
2614	lai_per_level_in, light_tran_to_level_in, gs_diffuco_output_in, &
2615	gstot_output_in)
2616
2617	USE netcdf
2618
2619	IMPLICIT NONE
2620
2621	!! 0.1 Input variables
2622
2623	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
2624
2625	CHARACTER(LEN=80) :: jamescdfname
2626
2627	REAL(r_std), DIMENSION (nvm, nlevels_tot), INTENT(in) :: gs_distribution_in
2628	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: gs_output_in
2629	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: gstot_output_in
2630	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: lai_per_level_in
2631	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: light_tran_to_level_in
2632	REAL(r_std),DIMENSION (nvm, nlevels_tot), INTENT(in) :: gs_diffuco_output_in
2633
2634
2635	LOGICAL :: james_netcdf_flag
2636
2637
2638	!! 0.2 Output variables
2639
2640	!! 0.3 Modified variables
2641
2642	!! 0.4 Local variables
2643
2644	INTEGER(i_std), SAVE :: gsd_id
2645	INTEGER(i_std), SAVE :: gso_id
2646	INTEGER(i_std), SAVE :: gsto_id
2647	INTEGER(i_std), SAVE :: lpl_id
2648	INTEGER(i_std), SAVE :: lttl_id
2649	INTEGER(i_std), SAVE :: gsdo_id
2650
2651	INTEGER(i_std) :: ii, jj
2652	INTEGER(i_std), SAVE :: james_step_count
2653
2654	INTEGER(i_std), SAVE :: iret, iret2, fid, fid2
2655	INTEGER(i_std) :: TimeDimID, jtid, ts_id, HeightDimID, i
2656	INTEGER(i_std) :: height_dim_id, height_dim_id2, time_dim_id
2657	INTEGER(i_std) :: pftno_dim_id
2658	INTEGER(i_std) :: dimids(3)
2659	INTEGER(i_std), SAVE :: tap_id
2660	INTEGER, DIMENSION(3) :: jstartnc, jcountnc
2661
2662
2663	jamescdfname = '/home/users/jryder/james.ryder/outputcdf.nc'
2664
2665	IF (james_netcdf_flag) THEN
2666
2667	james_netcdf_flag = .FALSE.
2668
2669	! james_step_count = 0
2670
2671	! WRITE(numout,*) 'james_netcdf_flag call here'
2672
2673	! âŠ define the parameters here
2674
2675	iret = nf90_create(jamescdfname, nf90_write, fid)
2676
2677
2678	! define the dimensions
2679
2680	iret = nf90_def_dim(fid, "pftno", nvm, pftno_dim_id)
2681	iret = nf90_def_dim(fid, "height", nlevels_tot, height_dim_id)
2682	iret = nf90_def_dim(fid, "time", nf90_unlimited, time_dim_id)
2683
2684	dimids = (/ pftno_dim_id, height_dim_id, time_dim_id /)
2685
2686
2687	! define the variables
2688	iret = nf90_def_var (fid, "gs_distribution_in", nf90_real4, dimids, gsd_id)
2689	iret = nf90_def_var (fid, "gs_output_in", nf90_real4, dimids, gso_id)
2690	iret = nf90_def_var (fid, "gstot_output_in", nf90_real4, dimids, gsto_id)
2691	iret = nf90_def_var (fid, "lai_per_level_in", nf90_real4, dimids, lpl_id)
2692	iret = nf90_def_var (fid, "light_tran_to_level_in", nf90_real4, dimids, lttl_id)
2693	iret = nf90_def_var (fid, "gs_diffuco_output_in", nf90_real4, dimids, gsdo_id)
2694
2695
2696	iret = nf90_enddef (fid)
2697
2698	END IF !(james_netcdf_flag)
2699
2700
2701	! open the netCDF file at each step
2702	iret = nf90_open(jamescdfname, nf90_Write, fid)
2703
2704	jstartnc = (/ 1, 1, 1 /) ! vector of integers specifying the index in the variable from
2705	! which the first of the data values will be read
2706	jcountnc = (/ nvm, nlevels_tot, 1 /) ! vector of integers specifying the number of indices
2707	! selected along each dimension
2708
2709	james_step_count = james_step_count + 1 ! this is the time step counter
2710
2711
2712	jstartnc(3) = james_step_count ! (2) refers to the time dimension at present (it
2713	! must always come last)
2714
2715
2716	iret = nf90_put_var(fid, gsd_id, gs_distribution_in, start = jstartnc, count = jcountnc)
2717	iret = nf90_put_var(fid, gso_id, gs_output_in, start = jstartnc, count = jcountnc)
2718	iret = nf90_put_var(fid, gsto_id, gstot_output_in, start = jstartnc, count = jcountnc)
2719	iret = nf90_put_var(fid, lpl_id, lai_per_level_in, start = jstartnc, count = jcountnc)
2720	iret = nf90_put_var(fid, lttl_id, light_tran_to_level_in, start = jstartnc, count = jcountnc)
2721	iret = nf90_put_var(fid, gsdo_id, gs_diffuco_output_in, start = jstartnc, count = jcountnc)
2722
2723
2724
2725	! close the netCDF file and write to disk
2726	iret = nf90_close(fid)
2727
2728
2729	END SUBROUTINE jnetcdf_trans_co2
2730
2731
2732	! -------------------------------------------------------------------------------------------------------------------------------------
2733
2734
2735
2736	END MODULE sechiba_hydrol_arch

Note: See TracBrowser for help on using the repository browser.

Download in other formats: