Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_sechiba/diffuco.f90 @ 7852

Last change on this file since 7852 was 7456, checked in by josefine.ghattas, 2 years ago
Add missing OMP THREADPRIVATE. Note this was not a significat error as the varaible is only used for printing.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 188.9 KB

Line
1	! ================================================================================================================================
2	! MODULE : diffuco
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.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 limiting coefficients, both aerodynamic
10	!! and hydrological, for the turbulent heat fluxes.
11	!!
12	!!\n DESCRIPTION: The aerodynamic resistance R_a is used to limit
13	!! the transport of fluxes from the surface layer of vegetation to the point in the atmosphere at which
14	!! interaction with the LMDZ atmospheric circulation model takes place. The aerodynamic resistance is
15	!! calculated either within the module r_aerod (if the surface drag coefficient is provided by the LMDZ, and
16	!! if the flag 'ldq_cdrag_from_gcm' is set to TRUE) or r_aero (if the surface drag coefficient must be calculated).\n
17	!!
18	!! Within ORCHIDEE, evapotranspiration is a function of the Evaporation Potential, but is modulated by a
19	!! series of resistances (canopy and aerodynamic) of the surface layer, here represented by beta.\n
20	!!
21	!! DESCRIPTION :
22	!! \latexonly
23	!! \input{diffuco_intro.tex}
24	!! \endlatexonly
25	!! \n
26	!!
27	!! This module calculates the beta for several different scenarios:
28	!! - diffuco_snow calculates the beta coefficient for sublimation by snow,
29	!! - diffuco_inter calculates the beta coefficient for interception loss by each type of vegetation,
30	!! - diffuco_bare calculates the beta coefficient for bare soil,
31	!! - diffuco_trans_co2 calculates the beta coefficient for transpiration for each type of vegetation, using Farqhar's formula
32	!! - chemistry_bvoc calculates the beta coefficient for emissions of biogenic compounds \n
33	!!
34	!! Finally, the module diffuco_comb computes the combined $\alpha$ and $\beta$ coefficients for use
35	!! elsewhere in the module. \n
36
37	!! RECENT CHANGE(S): Nathalie le 28 mars 2006 - sur proposition de Fred Hourdin, ajout
38	!! d'un potentiometre pour regler la resistance de la vegetation (rveg is now in pft_parameters)
39	!! October 2018: Removed diffuco_trans using Jarvis formula for calculation of beta coefficient
40	!!
41	!! REFERENCE(S) : None
42	!!
43	!! SVN :
44	!! $HeadURL$
45	!! $Date$
46	!! $Revision$
47	!! \n
48	!_ ================================================================================================================================
49
50	MODULE diffuco
51
52	! modules used :
53	USE constantes
54	USE constantes_soil
55	USE qsat_moisture
56	USE sechiba_io_p
57	USE ioipsl
58	USE pft_parameters
59	USE grid
60	USE time, ONLY : one_day, dt_sechiba, day_end
61	USE ioipsl_para
62	USE xios_orchidee
63	USE chemistry, ONLY : chemistry_initialize, chemistry_bvoc, chemistry_clear
64	USE function_library, ONLY: Arrhenius, Arrhenius_modified
65	IMPLICIT NONE
66
67	! public routines :
68	PRIVATE
69	PUBLIC :: diffuco_main, diffuco_initialize, diffuco_finalize, diffuco_clear
70
71	!
72	! variables used inside diffuco module : declaration and initialisation
73	!
74	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind module (m s^{-1})
75	!$OMP THREADPRIVATE(wind)
76	INTEGER(i_std), SAVE, PRIVATE :: printlev_loc !! Local level of text output for current module
77	!$OMP THREADPRIVATE(printlev_loc)
78
79
80	CONTAINS
81
82
83	!! =============================================================================================================================
84	!! SUBROUTINE: diffuco_initialize
85	!!
86	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
87	!!
88	!! DESCRIPTION: Allocate module variables, read from restart file or initialize with default values.
89	!! Call chemistry_initialize for initialization of variables needed for the calculations of BVOCs.
90	!!
91	!! RECENT CHANGE(S): None
92	!!
93	!! REFERENCE(S): None
94	!!
95	!! FLOWCHART: None
96	!! \n
97	!_ ==============================================================================================================================
98	SUBROUTINE diffuco_initialize (kjit, kjpindex, index, &
99	rest_id, lalo, neighbours, resolution, &
100	rstruct, q_cdrag)
101
102	!! 0. Variable and parameter declaration
103	!! 0.1 Input variables
104	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
105	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
106	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
107	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-)
108	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
109	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each
110	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
111
112	!! 0.2 Output variables
113	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! Structural resistance for the vegetation
114
115	!! 0.3 Modified variables
116	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
117
118	!! 0.4 Local variables
119	INTEGER :: ilai
120	INTEGER :: jv
121	INTEGER :: ier
122	CHARACTER(LEN=4) :: laistring
123	CHARACTER(LEN=80) :: var_name
124	!_ ================================================================================================================================
125
126	!! Initialize local printlev variable
127	printlev_loc=get_printlev('diffuco')
128
129
130	!! 1. Define flag ldq_cdrag_from_gcm. This flag determines if the cdrag should be taken from the GCM or be calculated.
131	!! The default value is true if the q_cdrag variables was already initialized. This is the case when coupled to the LMDZ.
132
133	!Config Key = CDRAG_FROM_GCM
134	!Config Desc = Keep cdrag coefficient from gcm.
135	!Config If = OK_SECHIBA
136	!Config Def = y
137	!Config Help = Set to .TRUE. if you want q_cdrag coming from GCM (if q_cdrag on initialization is non zero).
138	!Config Keep cdrag coefficient from gcm for latent and sensible heat fluxes.
139	!Config Units = [FLAG]
140	IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN
141	ldq_cdrag_from_gcm = .FALSE.
142	ELSE
143	ldq_cdrag_from_gcm = .TRUE.
144	ENDIF
145	CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm)
146	IF (printlev_loc>=2) WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm
147
148	!! 2. Allocate module variables
149	ALLOCATE (wind(kjpindex),stat=ier)
150	IF (ier /= 0) CALL ipslerr_p(3,'diffuco_initialize','Problem in allocate of variable wind','','')
151
152	!! 3. Read variables from restart file
153	IF (printlev>=3) WRITE (numout,*) 'Read DIFFUCO variables from restart file'
154
155	CALL ioconf_setatt_p('UNITS', 's/m')
156	CALL ioconf_setatt_p('LONG_NAME','Structural resistance')
157	CALL restget_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g)
158	IF ( ALL(rstruct(:,:) == val_exp) ) THEN
159	DO jv = 1, nvm
160	rstruct(:,jv) = rstruct_const(jv)
161	ENDDO
162	ENDIF
163
164	!! 4. Initialize chemistry module
165	IF (printlev>=3) WRITE(numout,*) "ok_bvoc:",ok_bvoc
166	IF ( ok_bvoc ) CALL chemistry_initialize(kjpindex, lalo, neighbours, resolution)
167
168	END SUBROUTINE diffuco_initialize
169
170
171
172	!! ================================================================================================================================
173	!! SUBROUTINE : diffuco_main
174	!!
175	!>\BRIEF The root subroutine for the module, which calls all other required
176	!! subroutines.
177	!!
178	!! DESCRIPTION :
179
180	!! This is the main subroutine for the module.
181	!! First it calculates the surface drag coefficient (via a call to diffuco_aero), using available parameters to determine
182	!! the stability of air in the surface layer by calculating the Richardson Nubmber. If a value for the
183	!! surface drag coefficient is passed down from the atmospheric model and and if the flag 'ldq_cdrag_from_gcm'
184	!! is set to TRUE, then the subroutine diffuco_aerod is called instead. This calculates the aerodynamic coefficient. \n
185	!!
186	!! Following this, an estimation of the saturated humidity at the surface is made (via a call
187	!! to qsatcalc in the module qsat_moisture). Following this the beta coefficients for sublimation (via
188	!! diffuco_snow), interception (diffuco_inter), bare soil (diffuco_bare), and transpiration (via
189	!! diffuco_trans_co2) are calculated in sequence. Finally
190	!! the alpha and beta coefficients are combined (diffuco_comb). \n
191	!!
192	!! The surface drag coefficient is calculated for use within the module enerbil. It is required to to
193	!! calculate the aerodynamic coefficient for every flux. \n
194	!!
195	!! The various beta coefficients are used within the module enerbil for modifying the rate of evaporation,
196	!! as appropriate for the surface. As explained in Chapter 2 of Guimberteau (2010), that module (enerbil)
197	!! calculates the rate of evaporation essentially according to the expression $E = /beta E_{pot}$, where
198	!! E is the total evaporation and $E_{pot}$ the evaporation potential. If $\beta = 1$, there would be
199	!! essentially no resistance to evaporation, whereas should $\beta = 0$, there would be no evaporation and
200	!! the surface layer would be subject to some very stong hydrological stress. \n
201	!!
202	!! The following processes are calculated:
203	!! - call diffuco_aero for aerodynamic transfer coeficient
204	!! - call diffuco_snow for partial beta coefficient: sublimation
205	!! - call diffuco_inter for partial beta coefficient: interception for each type of vegetation
206	!! - call diffuco_bare for partial beta coefficient: bare soil
207	!! - call diffuco_trans_co2 for partial beta coefficient: transpiration for each type of vegetation, using Farqhar's formula
208	!! - call diffuco_comb for alpha and beta coefficient
209	!! - call chemistry_bvoc for alpha and beta coefficients for biogenic emissions
210	!!
211	!! RECENT CHANGE(S): None
212	!!
213	!! MAIN OUTPUT VARIABLE(S): humrel, q_cdrag, vbeta, vbeta1, vbeta4,
214	!! vbeta2, vbeta3, rveget, cimean
215	!!
216	!! REFERENCE(S) :
217	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
218	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
219	!! Circulation Model. Journal of Climate, 6, pp.248-273.
220	!! - de Rosnay, P, 1999. ReprÃ©sentation des interactions sol-plante-atmosphÃšre dans le modÃšle de circulation gÃ©nÃ©rale
221	!! du LMD, 1999. PhD Thesis, UniversitÃ© Paris 6, available (25/01/12):
222	!! http://www.ecmwf.int/staff/patricia_de_rosnay/publications.html#8
223	!! - Ducharne, A, 1997. Le cycle de l'eau: modÃ©lisation de l'hydrologie continentale, Ã©tude de ses interactions avec
224	!! le climat, PhD Thesis, UniversitÃ© Paris 6
225	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
226	!! sur le cycle de l'eau, PhD Thesis, available (25/01/12):
227	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
228	!! - LathiÃšre, J, 2005. Evolution des Ã©missions de composÃ©s organiques et azotÃ©s par la biosphÃšre continentale dans le
229	!! modÃšle LMDz-INCA-ORCHIDEE, UniversitÃ© Paris 6
230	!!
231	!! FLOWCHART :
232	!! \latexonly
233	!! \includegraphics[scale=0.5]{diffuco_main_flowchart.png}
234	!! \endlatexonly
235	!! \n
236	!_ ================================================================================================================================
237
238	SUBROUTINE diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
239	& zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, q_cdrag, qsurf, qair, pb, &
240	& evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, &
241	& frac_nobio, snow_nobio, totfrac_nobio, &
242	& swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, &
243	& qsintveg, qsintmax, assim_param, &
244	& vbeta , vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, cresist, &
245	& lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
246	& hist_id, hist2_id, vbeta2sum, vbeta3sum, Light_Abs_Tot, &
247	& Light_Tran_Tot, lai_per_level, vbeta23, leaf_ci, &
248	& gs_distribution, gs_diffuco_output, gstot_component, gstot_frac,&
249	& warnings, u_speed, profile_vbeta3, profile_rveget, &
250	& delta_c13_assim, leaf_ci_out, info_limitphoto, JJ_out, assimi_lev)
251
252	!! 0. Variable and parameter declaration
253
254	!! 0.1 Input variables
255
256	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
257	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
258	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-)
259	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (-)
260	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
261	INTEGER(i_std),DIMENSION (kjpindex*(nlevels_tot+1)), INTENT (in) :: indexlai !! Indeces of the points on the 3D map
262	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-)
263	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
264	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
265	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
266	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m)
267	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m)
268	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness (m)
269	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
270	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
271	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
272	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
273	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface air specific humidity (kg kg^{-1})
274	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level air specific humidity (kg kg^{-1})
275	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
276	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! Fraction of floodplains
277	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! Reservoir in floodplains (estimation to avoid over-evaporation)
278	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure (hPa)
279	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! Limit to the bare soil evaporation when the
280	!! 11-layer hydrology is used (-)
281	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns!! Limit to the bare soil evaporation when the
282	!! 11-layer hydrology is used (-)
283	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1})
284	!! NdN This variable does not seem to be used at
285	!! all in diffuco
286	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation
287	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... (-)
288	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... (kg m^{-2})
289	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... (-)
290	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
291	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux (W m^{-2})
292	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless)
293	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
294	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (-)
295	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
296	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! PFT leaf area index (m^{2} m^{-2})
297	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
298	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
299	!! (kg m^{-2})
300	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! referenced vcmax, nue, leaf N for photosynthesis
301	!! Nitrogen use Efficiency with impact of leaf age (umol CO2 (gN)-1 s-1)
302	!! Total canopy (eg leaf) N (gN m-2[ground])
303	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
304	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each
305	!! grid point (1=N, 2=E, 3=S, 4=W)
306	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
307	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ptnlev1 !! 1st level of soil temperature (Kelvin)
308	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation expressed in mm/tstep
309	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (in) :: frac_age !! Age efficiency for isoprene emissions (from STOMATE)
310	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
311	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: Light_Abs_Tot !! Absorbed radiation per level for photosynthesis
312	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: Light_Tran_Tot !! Transmitted radiation per level for photosynthesis
313	REAL(r_std), DIMENSION(:,:,:), INTENT(IN) :: lai_per_level !! This is the LAI per vertical level
314	!! @tex $(m^{2} m^{-2})$
315	REAL(r_std), DIMENSION (:), INTENT (in) :: u_speed !! wind speed at specified canopy level (m/s)
316	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
317	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
318
319
320	!! 0.2 Output variables
321
322	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
323	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (-)
324	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
325	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains
326	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! Mean stomatal conductance to CO2 (mol m-2 s-1)
327	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
328	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration (-)
329	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for potential transpiration
330	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Stomatal resistance for the whole canopy (s m^{-1})
331	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: rstruct !! Structural resistance for the vegetation
332	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! Mean leaf Ci (ppm)
333	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area)
334	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cresist !! coefficient for resistances (??)
335
336	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will
337	!! transpire once intercepted water has evaporated (-)
338	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: leaf_ci !! intercellular CO2 concentration (ppm)
339	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gs_distribution
340	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gs_diffuco_output
341	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gstot_component
342	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: gstot_frac
343	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: profile_vbeta3
344	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: profile_rveget
345
346	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: delta_c13_assim !! C13 concentration in delta notation
347	!! @tex $ permille $ @endtex (per thousand)
348	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: leaf_ci_out !! Ci
349	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot+1),INTENT(out) :: info_limitphoto
350	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot),INTENT(out) :: JJ_out
351	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot),INTENT(out) :: assimi_lev !! assimilation (at all LAI levels)
352	!! ($\mumolm^{-2}s^{-1}$)
353
354	!! 0.3 Modified variables
355
356	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
357	REAL(r_std),DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) &
358	:: warnings !! A warning counter
359	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag coefficient (-)
360
361
362	!! 0.4 Local variables
363	REAL(r_std),DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum
364	REAL(r_std),DIMENSION(kjpindex) :: raero !! Aerodynamic resistance (s m^{-1})
365	INTEGER(i_std) :: ilai
366	CHARACTER(LEN=4) :: laistring
367	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
368	REAL(r_std),DIMENSION(kjpindex) :: qsatt !! Surface saturated humidity (kg kg^{-1})
369	REAL(r_std),DIMENSION (kjpindex,nvm) :: cim !! Intercellular CO2 over nlevels_tot
370
371	! Isotope_c13
372	INTEGER(i_std), DIMENSION(kjpindex) :: index_assi !! indices (unitless)
373	INTEGER(i_std) :: nia,inia,nina,iainia !! counter/indices (unitless)
374	INTEGER(i_std) :: ji,jl,jv !! Indices (-)
375	REAL(r_std), DIMENSION(nlevels_tot+1) :: laitab !! LAI per layer (m^2.m^{-2})
376	REAL(r_std), DIMENSION(kjpindex,nvm) :: assimtot !! total assimilation
377	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
378
379	!_ ================================================================================================================================
380
381
382	gs_distribution = 0.0d0
383	wind(:) = SQRT (u(:)u(:) + v(:)v(:))
384
385	!! 1. Calculate the different coefficients
386
387	IF (.NOT.ldq_cdrag_from_gcm) THEN
388	! Case 1a)
389	CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
390	qsurf, qair, snow, q_cdrag)
391	ENDIF
392
393	! Case 1b)
394	CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
395
396	!! 2. Make an estimation of the saturated humidity at the surface
397
398	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
399
400	!! 3. Calculate the beta coefficient for sublimation
401
402	CALL diffuco_snow (kjpindex, qair, qsatt, rau, u, v, q_cdrag, &
403	snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, &
404	vbeta1)
405
406	CALL diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
407	& flood_frac, flood_res, vbeta5)
408
409	!! 4. Calculate the beta coefficient for interception
410
411	CALL diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
412	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
413
414	!! 5. Calculate the beta coefficient for transpiration
415	CALL diffuco_trans_co2 (kjpindex, swdown, pb, qair, temp_air, &
416	temp_growth, rau, q_cdrag, humrel, &
417	assim_param, ccanopy, lai, &
418	veget_max, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, &
419	rstruct, cimean, gsmean, gpp, cresist, vbeta23, Light_Abs_Tot,&
420	Light_Tran_Tot, lai_per_level, gs_distribution, gs_diffuco_output, &
421	gstot_component, gstot_frac, veget, warnings, u_speed, profile_vbeta3, &
422	profile_rveget, index_assi, leaf_ci, assimtot, assimi_lev, hist_id, &
423	indexveg, indexlai, index, kjit, cim, info_limitphoto, JJ_out)
424
425	! calculate C13 isotopic signature
426	IF (ok_c13) THEN
427	CALL isotope_c13(kjpindex, index_assi, leaf_ci, ccanopy, lai_per_level, &
428	assimtot, delta_c13_assim, leaf_ci_out, nvm, &
429	nlevels_tot, assimi_lev)
430	ENDIF
431
432	!
433	!biogenic emissions
434	!
435
436	IF ( ok_bvoc ) THEN
437	CALL chemistry_bvoc (kjpindex, swdown, coszang, temp_air, &
438	temp_sol, ptnlev1, precip_rain, humrel, veget_max, &
439	lai, frac_age, lalo, ccanopy, cim, wind, snow, &
440	veget, hist_id, hist2_id, kjit, index, &
441	indexlai, indexveg)
442	ENDIF
443	!
444	! combination of coefficient : alpha and beta coefficient
445	! beta coefficient for bare soil
446	!
447
448	CALL diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)
449
450	!! 6. Combine the alpha and beta coefficients
451
452	CALL diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
453	snow, veget, lai, &
454	tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, &
455	evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax, vbeta2sum, vbeta3sum)
456
457	CALL xios_orchidee_send_field("q_cdrag",q_cdrag)
458	CALL xios_orchidee_send_field("raero",raero)
459	CALL xios_orchidee_send_field("wind",wind)
460	CALL xios_orchidee_send_field("qsatt",qsatt)
461	CALL xios_orchidee_send_field("coszang",coszang)
462	CALL xios_orchidee_send_field('cim', cim)
463
464	IF ( .NOT. almaoutput ) THEN
465	CALL histwrite_p(hist_id, 'raero', kjit, raero, kjpindex, index)
466	CALL histwrite_p(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
467	CALL histwrite_p(hist_id, 'Wind', kjit, wind, kjpindex, index)
468	CALL histwrite_p(hist_id, 'qsatt', kjit, qsatt, kjpindex, index)
469	CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)
470
471	IF ( hist2_id > 0 ) THEN
472	CALL histwrite_p(hist2_id, 'raero', kjit, raero, kjpindex, index)
473	CALL histwrite_p(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
474	CALL histwrite_p(hist2_id, 'Wind', kjit, wind, kjpindex, index)
475	CALL histwrite_p(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index)
476	ENDIF
477	ELSE
478	CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)
479	ENDIF
480
481	IF (printlev>=3) WRITE (numout,*) ' diffuco_main done '
482
483	END SUBROUTINE diffuco_main
484
485	!! =============================================================================================================================
486	!! SUBROUTINE: diffuco_finalize
487	!!
488	!>\BRIEF Write to restart file
489	!!
490	!! DESCRIPTION: This subroutine writes the module variables and variables calculated in diffuco
491	!! to restart file
492	!!
493	!! RECENT CHANGE(S): None
494	!! REFERENCE(S): None
495	!! FLOWCHART: None
496	!! \n
497	!_ ==============================================================================================================================
498	SUBROUTINE diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
499
500	!! 0. Variable and parameter declaration
501	!! 0.1 Input variables
502	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
503	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
504	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-)
505	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! Structural resistance for the vegetation
506
507	!! 0.4 Local variables
508	INTEGER :: ilai
509	CHARACTER(LEN=4) :: laistring
510	CHARACTER(LEN=80) :: var_name
511
512	!_ ================================================================================================================================
513
514	!! 1. Prepare the restart file for the next simulation
515	IF (printlev>=3) WRITE (numout,*) 'Complete restart file with DIFFUCO variables '
516
517	CALL restput_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g)
518
519	END SUBROUTINE diffuco_finalize
520
521
522	!! ================================================================================================================================
523	!! SUBROUTINE : diffuco_clear
524	!!
525	!>\BRIEF Housekeeping module to deallocate the variables
526	!! rstruct and raero
527	!!
528	!! DESCRIPTION : Housekeeping module to deallocate the variables
529	!! rstruct and raero
530	!!
531	!! RECENT CHANGE(S) : None
532	!!
533	!! MAIN OUTPUT VARIABLE(S) : None
534	!!
535	!! REFERENCE(S) : None
536	!!
537	!! FLOWCHART : None
538	!! \n
539	!_ ================================================================================================================================
540
541	SUBROUTINE diffuco_clear()
542
543	! Deallocate and reset variables in chemistry module
544	CALL chemistry_clear
545
546	END SUBROUTINE diffuco_clear
547
548
549	!! ================================================================================================================================
550	!! SUBROUTINE : diffuco_aero
551	!!
552	!>\BRIEF This module first calculates the surface drag
553	!! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled
554	!! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE
555	!!
556	!! DESCRIPTION : Computes the surface drag coefficient, for cases
557	!! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the
558	!! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric
559	!! stability in the surface layer. The formulation used to find this surface drag coefficient is
560	!! dependent on the stability determined. \n
561	!!
562	!! Designation of wind speed
563	!! \latexonly
564	!! \input{diffucoaero1.tex}
565	!! \endlatexonly
566	!!
567	!! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson
568	!! eqn.4.47, 2005). (required for calculation of the Richardson Number)
569	!! \latexonly
570	!! \input{diffucoaero2.tex}
571	!! \endlatexonly
572	!!
573	!! \latexonly
574	!! \input{diffucoaero3.tex}
575	!! \endlatexonly
576	!!
577	!! Calculation of the virtual air temperature at the surface (required for calculation
578	!! of the Richardson Number)
579	!! \latexonly
580	!! \input{diffucoaero4.tex}
581	!! \endlatexonly
582	!!
583	!! Calculation of the virtual surface temperature (required for calculation of th
584	!! Richardson Number)
585	!! \latexonly
586	!! \input{diffucoaero5.tex}
587	!! \endlatexonly
588	!!
589	!! Calculation of the squared wind shear (required for calculation of the Richardson
590	!! Number)
591	!! \latexonly
592	!! \input{diffucoaero6.tex}
593	!! \endlatexonly
594	!!
595	!! Calculation of the Richardson Number. The Richardson Number is defined as the ratio
596	!! of potential to kinetic energy, or, in the context of atmospheric science, of the
597	!! generation of energy by wind shear against consumption
598	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
599	!! becomes turbulent and vise versa). It is approximated using the expression below:
600	!! \latexonly
601	!! \input{diffucoaero7.tex}
602	!! \endlatexonly
603	!!
604	!! The Richardson Number hence calculated is subject to a minimum value:
605	!! \latexonly
606	!! \input{diffucoaero8.tex}
607	!! \endlatexonly
608	!!
609	!! Computing the drag coefficient. We add the add the height of the vegetation to the
610	!! level height to take into account that the level 'seen' by the vegetation is actually
611	!! the top of the vegetation. Then we we can subtract the displacement height.
612	!! \latexonly
613	!! \input{diffucoaero9.tex}
614	!! \endlatexonly
615	!!
616	!! For the stable case (i.e $R_i$ $\geq$ 0)
617	!! \latexonly
618	!! \input{diffucoaero10.tex}
619	!! \endlatexonly
620	!!
621	!! \latexonly
622	!! \input{diffucoaero11.tex}
623	!! \endlatexonly
624	!!
625	!! For the unstable case (i.e. $R_i$ < 0)
626	!! \latexonly
627	!! \input{diffucoaero12.tex}
628	!! \endlatexonly
629	!!
630	!! \latexonly
631	!! \input{diffucoaero13.tex}
632	!! \endlatexonly
633	!!
634	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
635	!! To prevent this, a minimum limit to the drag coefficient is defined as:
636	!!
637	!! \latexonly
638	!! \input{diffucoaero14.tex}
639	!! \endlatexonly
640	!!
641	!! RECENT CHANGE(S): None
642	!!
643	!! MAIN OUTPUT VARIABLE(S): q_cdrag
644	!!
645	!! REFERENCE(S) :
646	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
647	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
648	!! Circulation Model. Journal of Climate, 6, pp.248-273
649	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
650	!! sur le cycle de l'eau, PhD Thesis, available from:
651	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
652	!! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge
653	!! University Press, ISBN 0-521-54865-9
654	!!
655	!! FLOWCHART :
656	!! \latexonly
657	!! \includegraphics[scale=0.5]{diffuco_aero_flowchart.png}
658	!! \endlatexonly
659	!! \n
660	!_ ================================================================================================================================
661
662	SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
663	qsurf, qair, snow, q_cdrag)
664
665	!! 0. Variable and parameter declaration
666
667	!! 0.1 Input variables
668
669	INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size
670	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
671	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
672	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first atmospheric layer (m)
673	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m)
674	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m)
675	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height (m)
676	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Ground temperature (K)
677	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
678	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific air humidity (kg kg^{-1})
679	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
680	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
681
682	!! 0.2 Output variables
683
684	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_cdrag !! Surface drag coefficient (-)
685
686	!! 0.3 Modified variables
687
688	!! 0.4 Local variables
689
690	INTEGER(i_std) :: ji, jv
691	REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2
692	REAL(r_std) :: zri, cd_neut, zscf, cd_tmp
693	!_ ================================================================================================================================
694
695	!! 1. Initialisation
696
697	! test if we have to work with q_cdrag or to calcul it
698	DO ji=1,kjpindex
699
700	!! 1a).1 Designation of wind speed
701	!! \latexonly
702	!! \input{diffucoaero1.tex}
703	!! \endlatexonly
704	speed = wind(ji)
705
706	!! 1a).2 Calculation of geopotentiel
707	!! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required
708	!! for calculation of the Richardson Number)
709	!! \latexonly
710	!! \input{diffucoaero2.tex}
711	!! \endlatexonly
712	zg = zlev(ji) * cte_grav
713
714	!! \latexonly
715	!! \input{diffucoaero3.tex}
716	!! \endlatexonly
717	zdphi = zg/cp_air
718
719	!! 1a).3 Calculation of the virtual air temperature at the surface
720	!! required for calculation of the Richardson Number
721	!! \latexonly
722	!! \input{diffucoaero4.tex}
723	!! \endlatexonly
724	ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji))
725
726	!! 1a).4 Calculation of the virtual surface temperature
727	!! required for calculation of the Richardson Number
728	!! \latexonly
729	!! \input{diffucoaero5.tex}
730	!! \endlatexonly
731	ztvs = temp_sol(ji) * (un + retv * qsurf(ji))
732
733	!! 1a).5 Calculation of the squared wind shear
734	!! required for calculation of the Richardson Number
735	!! \latexonly
736	!! \input{diffucoaero6.tex}
737	!! \endlatexonly
738	zdu2 = MAX(cepdu2,speed**2)
739
740	!! 1a).6 Calculation of the Richardson Number
741	!! The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the
742	!! context of atmospheric science, of the generation of energy by wind shear against consumption
743	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
744	!! becomes turbulent and vise versa).\n
745	!! It is approximated using the expression below:
746	!! \latexonly
747	!! \input{diffucoaero7.tex}
748	!! \endlatexonly
749	zri = zg * (ztvd - ztvs) / (zdu2 * ztvd)
750
751	!! The Richardson Number hence calculated is subject to a minimum value:
752	!! \latexonly
753	!! \input{diffucoaero8.tex}
754	!! \endlatexonly
755	zri = MAX(MIN(zri,5.),-5.)
756
757	!! 1a).7 Computing the drag coefficient
758	!! We add the add the height of the vegetation to the level height to take into account
759	!! that the level 'seen' by the vegetation is actually the top of the vegetation. Then we
760	!! we can subtract the displacement height.
761	!! \latexonly
762	!! \input{diffucoaero9.tex}
763	!! \endlatexonly
764
765	!! 7.0 Snow smoothering
766	!! Snow induces low levels of turbulence.
767	!! Sensible heat fluxes can therefore be reduced of ~1/3. Pomeroy et al., 1998
768	cd_neut = ct_karman ** 2. / ( LOG( (zlev(ji) + roughheight(ji)) / z0m(ji) ) * LOG( (zlev(ji) + roughheight(ji)) / z0h(ji) ) )
769
770	!! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0)
771	IF (zri .GE. zero) THEN
772
773	!! \latexonly
774	!! \input{diffucoaero10.tex}
775	!! \endlatexonly
776	zscf = SQRT(un + cd * ABS(zri))
777
778	!! \latexonly
779	!! \input{diffucoaero11.tex}
780	!! \endlatexonly
781	cd_tmp=cd_neut/(un + trois * cb * zri * zscf)
782	ELSE
783
784	!! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0)
785	!! \latexonly
786	!! \input{diffucoaero12.tex}
787	!! \endlatexonly
788	zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * &
789	& ((zlev(ji) + roughheight(ji)) / z0m(ji))))
790
791	!! \latexonly
792	!! \input{diffucoaero13.tex}
793	!! \endlatexonly
794	cd_tmp=cd_neut * (un - trois * cb * zri * zscf)
795	ENDIF
796
797	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
798	!! To prevent this, a minimum limit to the drag coefficient is defined as:
799
800	!! \latexonly
801	!! \input{diffucoaero14.tex}
802	!! \endlatexonly
803	!!
804	q_cdrag(ji) = MAX(cd_tmp, min_qc/MAX(speed,min_wind))
805
806	! In some situations it might be useful to give an upper limit on the cdrag as well.
807	! The line here should then be uncommented.
808	!q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind))
809
810	END DO
811
812	IF (printlev_loc>=3) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done '
813
814	END SUBROUTINE diffuco_aero
815
816
817	!! ================================================================================================================================
818	!! SUBROUTINE : diffuco_snow
819	!!
820	!>\BRIEF This subroutine computes the beta coefficient for snow sublimation.
821	!!
822	!! DESCRIPTION : This routine computes beta coefficient for snow sublimation, which
823	!! integrates the snow on both vegetation and other surface types (e.g. ice, lakes,
824	!! cities etc.) \n
825	!!
826	!! A critical depth of snow (snowcri) is defined to calculate the fraction of each grid-cell
827	!! that is covered with snow (snow/snowcri) while the remaining part is snow-free.
828	!! We also carry out a first calculation of sublimation (subtest) to lower down the beta
829	!! coefficient if necessary (if subtest > snow). This is a predictor-corrector test.
830	!!
831	!! RECENT CHANGE(S): None
832	!!
833	!! MAIN OUTPUT VARIABLE(S): ::vbeta1
834	!!
835	!! REFERENCE(S) :
836	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
837	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
838	!! Circulation Model. Journal of Climate, 6, pp. 248-273
839	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
840	!! sur le cycle de l'eau, PhD Thesis, available from:
841	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
842	!!
843	!! FLOWCHART : None
844	!! \n
845	!_ ================================================================================================================================
846
847	SUBROUTINE diffuco_snow (kjpindex, qair, qsatt, rau, u, v,q_cdrag, &
848	& snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, &
849	vbeta1)
850
851	!! 0. Variable and parameter declaration
852
853	!! 0.1 Input variables
854
855	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
856	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
857	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
858	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
859	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
860	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
861	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
862	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg m^{-2})
863	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. (-)
864	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice, lakes, cities etc. (-)
865	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice, lakes, cities etc. (-)
866	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
867	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
868
869	!! 0.2 Output variables
870
871	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (dimensionless ratio)
872
873	!! 0.3 Modified variables
874
875	!! 0.4 Local variables
876
877	REAL(r_std) :: subtest !! Sublimation for test (kg m^{-2})
878	REAL(r_std) :: zrapp !! Modified factor (ratio)
879	REAL(r_std) :: speed !! Wind speed (m s^{-1})
880	REAL(r_std) :: vbeta1_add !! Beta for sublimation (ratio)
881	INTEGER(i_std) :: ji, jv !! Indices (-)
882	!_ ================================================================================================================================
883
884	!! 1. Calculate beta coefficient for snow sublimation on the vegetation\n
885
886	DO ji=1,kjpindex ! Loop over # pixels - domain size
887
888	! Fraction of mesh that can sublimate snow
889	vbeta1(ji) = (un - totfrac_nobio(ji)) * frac_snow_veg(ji)
890
891	! Limitation of sublimation in case of snow amounts smaller
892	! than the atmospheric demand.
893	speed = MAX(min_wind, wind(ji))
894
895	subtest = dt_sechiba * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * &
896	& ( qsatt(ji) - qair(ji) )
897
898	IF ( subtest .GT. min_sechiba ) THEN
899	zrapp = snow(ji) / subtest
900	IF ( zrapp .LT. un ) THEN
901	vbeta1(ji) = vbeta1(ji) * zrapp
902	ENDIF
903	ENDIF
904
905	END DO ! Loop over # pixels - domain size
906
907	!! 2. Add the beta coefficients calculated from other surfaces types (snow on ice,lakes, cities...)
908
909	DO jv = 1, nnobio ! Loop over # other surface types
910	!!$ !
911	!!$ IF ( jv .EQ. iice ) THEN
912	!!$ !
913	!!$ ! Land ice is of course a particular case
914	!!$ !
915	!!$ DO ji=1,kjpindex
916	!!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv)
917	!!$ ENDDO
918	!!$ !
919	!!$ ELSE
920	!
921	DO ji=1,kjpindex ! Loop over # pixels - domain size
922
923	vbeta1_add = frac_nobio(ji,jv) * frac_snow_nobio(ji, jv)
924
925	! Limitation of sublimation in case of snow amounts smaller than
926	! the atmospheric demand.
927	speed = MAX(min_wind, wind(ji))
928
929	!! Limitation of sublimation by the snow accumulated on the ground
930	!! A first approximation is obtained with the old values of
931	!! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc)
932	subtest = dt_sechiba * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * &
933	& ( qsatt(ji) - qair(ji) )
934
935	IF ( subtest .GT. min_sechiba ) THEN
936	zrapp = snow_nobio(ji,jv) / subtest
937	IF ( zrapp .LT. un ) THEN
938	vbeta1_add = vbeta1_add * zrapp
939	ENDIF
940	ENDIF
941
942	vbeta1(ji) = vbeta1(ji) + vbeta1_add
943
944	ENDDO ! Loop over # pixels - domain size
945
946	!!$ ENDIF
947
948	ENDDO ! Loop over # other surface types
949
950	IF (printlev>=3) WRITE (numout,*) ' diffuco_snow done '
951
952	END SUBROUTINE diffuco_snow
953
954
955	!! ================================================================================================================================
956	!! SUBROUTINE : diffuco_flood
957	!!
958	!>\BRIEF This routine computes partial beta coefficient : floodplains
959	!!
960	!! DESCRIPTION :
961	!!
962	!! RECENT CHANGE(S): None
963	!!
964	!! MAIN OUTPUT VARIABLE(S) : vbeta5
965	!!
966	!! REFERENCE(S) : None
967	!!
968	!! FLOWCHART : None
969	!! \n
970	!_ ================================================================================================================================
971
972	SUBROUTINE diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
973	& flood_frac, flood_res, vbeta5)
974
975	! interface description
976	! input scalar
977	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
978	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
979	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity
980	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
981	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
982	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
983	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
984	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! water mass in flood reservoir
985	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! fraction of floodplains
986	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Potential evaporation
987	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr!! Potential evaporation2
988	! output fields
989	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains
990
991	! local declaration
992	REAL(r_std) :: subtest, zrapp, speed
993	INTEGER(i_std) :: ji, jv
994
995	!_ ================================================================================================================================
996	!
997	! beta coefficient for sublimation for floodplains
998	!
999	DO ji=1,kjpindex
1000	!
1001	IF (evapot(ji) .GT. min_sechiba) THEN
1002	vbeta5(ji) = flood_frac(ji) *evapot_corr(ji)/evapot(ji)
1003	ELSE
1004	vbeta5(ji) = flood_frac(ji)
1005	ENDIF
1006	!
1007	! -- Limitation of evaporation in case of water amounts smaller than
1008	! the atmospheric demand.
1009
1010	!
1011	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
1012	!
1013	subtest = dt_sechiba * vbeta5(ji) * speed * q_cdrag(ji) * rau(ji) * &
1014	& ( qsatt(ji) - qair(ji) )
1015	!
1016	IF ( subtest .GT. min_sechiba ) THEN
1017	zrapp = flood_res(ji) / subtest
1018	IF ( zrapp .LT. un ) THEN
1019	vbeta5(ji) = vbeta5(ji) * zrapp
1020	ENDIF
1021	ENDIF
1022	!
1023	END DO
1024
1025	IF (printlev>=3) WRITE (numout,*) ' diffuco_flood done '
1026
1027	END SUBROUTINE diffuco_flood
1028
1029
1030	!! ================================================================================================================================
1031	!! SUBROUTINE : diffuco_inter
1032	!!
1033	!>\BRIEF This routine computes the partial beta coefficient
1034	!! for the interception for each type of vegetation
1035	!!
1036	!! DESCRIPTION : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax).
1037	!! The former is submitted to transpiration only (vbeta3 coefficient, calculated in
1038	!! diffuco_trans_co2), while the latter is first submitted to interception loss
1039	!! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated
1040	!! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test,
1041	!! as for snow sublimation. \n
1042	!!
1043	!! \latexonly
1044	!! \input{diffucointer1.tex}
1045	!! \endlatexonly
1046	!! Calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1047	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1048	!!
1049	!! \latexonly
1050	!! \input{diffucointer2.tex}
1051	!! \endlatexonly
1052	!!
1053	!! Calculation of $\beta_3$, the canopy transpiration resistance
1054	!! \latexonly
1055	!! \input{diffucointer3.tex}
1056	!! \endlatexonly
1057	!!
1058	!! We here determine the limitation of interception loss by the water stored on the leaf.
1059	!! A first approximation of interception loss is obtained using the old values of
1060	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1061	!! \latexonly
1062	!! \input{diffucointer4.tex}
1063	!! \endlatexonly
1064	!!
1065	!! \latexonly
1066	!! \input{diffucointer5.tex}
1067	!! \endlatexonly
1068	!!
1069	!! \latexonly
1070	!! \input{diffucointer6.tex}
1071	!! \endlatexonly
1072	!!
1073	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1074	!! 'zqsvegrap'.
1075	!! \latexonly
1076	!! \input{diffucointer7.tex}
1077	!! \endlatexonly
1078	!!
1079	!! RECENT CHANGE(S): None
1080	!!
1081	!! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23
1082	!!
1083	!! REFERENCE(S) :
1084	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1085	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1086	!! Circulation Model. Journal of Climate, 6, pp. 248-273
1087	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1088	!! sur le cycle de l'eau, PhD Thesis, available from:
1089	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1090	!! - Perrier, A, 1975. Etude physique de l'Ã©vaporation dans les conditions naturelles. Annales
1091	!! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243
1092	!!
1093	!! FLOWCHART : None
1094	!! \n
1095	!_ ================================================================================================================================
1096
1097	SUBROUTINE diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
1098	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
1099
1100	!! 0 Variable and parameter declaration
1101
1102	!! 0.1 Input variables
1103
1104	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1105	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
1106	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
1107	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
1108	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
1109	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
1110	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1111	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
1112	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type (fraction)
1113	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
1114	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
1115	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! architectural resistance (s m^{-1})
1116
1117	!! 0.2 Output variables
1118
1119	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
1120	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will
1121	!! transpire (-)
1122
1123	!! 0.4 Local variables
1124
1125	INTEGER(i_std) :: ji, jv !! (-), (-)
1126	REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed !!
1127	!_ ================================================================================================================================
1128
1129	!! 1. Initialize
1130
1131	vbeta2(:,:) = zero
1132	vbeta23(:,:) = zero
1133
1134	!! 2. The beta coefficient for interception by vegetation.
1135
1136	DO jv = 2,nvm
1137
1138	DO ji=1,kjpindex
1139
1140	IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN
1141
1142	zqsvegrap = zero
1143	IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN
1144
1145	!! \latexonly
1146	!! \input{diffucointer1.tex}
1147	!! \endlatexonly
1148	!!
1149	!! We calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1150	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1151	!!
1152	zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
1153	END IF
1154
1155	!! \latexonly
1156	!! \input{diffucointer2.tex}
1157	!! \endlatexonly
1158	speed = MAX(min_wind, wind(ji))
1159
1160	!! Calculation of $\beta_3$, the canopy transpiration resistance
1161	!! \latexonly
1162	!! \input{diffucointer3.tex}
1163	!! \endlatexonly
1164	vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv)))
1165
1166	!! We here determine the limitation of interception loss by the water stored on the leaf.
1167	!! A first approximation of interception loss is obtained using the old values of
1168	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1169	!! \latexonly
1170	!! \input{diffucointer4.tex}
1171	!! \endlatexonly
1172	ziltest = dt_sechiba * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * &
1173	& ( qsatt(ji) - qair(ji) )
1174
1175	IF ( ziltest .GT. min_sechiba ) THEN
1176
1177	!! \latexonly
1178	!! \input{diffucointer5.tex}
1179	!! \endlatexonly
1180	zrapp = qsintveg(ji,jv) / ziltest
1181	IF ( zrapp .LT. un ) THEN
1182
1183	!! \latexonly
1184	!! \input{diffucointer6.tex}
1185	!! \endlatexonly
1186	!!
1187	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1188	!! 'zqsvegrap'.
1189	IF ( humrel(ji,jv) >= min_sechiba ) THEN
1190	vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, zero)
1191	ELSE
1192	! We don't want transpiration when the soil cannot deliver it
1193	vbeta23(ji,jv) = zero
1194	ENDIF
1195
1196	!! \latexonly
1197	!! \input{diffucointer7.tex}
1198	!! \endlatexonly
1199	vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp
1200	ENDIF
1201	ENDIF
1202	END IF
1203	! ! Autre formulation possible pour l'evaporation permettant une transpiration sur tout le feuillage
1204	! !commenter si formulation Nathalie sinon Tristan
1205	! speed = MAX(min_wind, wind(ji))
1206	!
1207	! vbeta23(ji,jv) = MAX(zero, veget(ji,jv) * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) - vbeta2(ji,jv))
1208
1209	END DO
1210
1211	END DO
1212
1213	IF (printlev>=3) WRITE (numout,*) ' diffuco_inter done '
1214
1215	END SUBROUTINE diffuco_inter
1216
1217
1218	!! ==============================================================================================================================
1219	!! SUBROUTINE : diffuco_bare
1220	!!
1221	!>\BRIEF This routine computes the partial beta coefficient corresponding to
1222	!! bare soil
1223	!!
1224	!! DESCRIPTION : Bare soil evaporation is submitted to a maximum possible flow (evap_bare_lim)
1225	!!
1226	!! Calculation of wind speed
1227	!! \latexonly
1228	!! \input{diffucobare1.tex}
1229	!! \endlatexonly
1230	!!
1231	!! The calculation of $\beta_4$
1232	!! \latexonly
1233	!! \input{diffucobare2.tex}
1234	!! \endlatexonly
1235	!!
1236	!! RECENT CHANGE(S): None
1237	!!
1238	!! MAIN OUTPUT VARIABLE(S): ::vbeta4
1239	!!
1240	!! REFERENCE(S) :
1241	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1242	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1243	!! Circulation Model. Journal of Climate, 6, pp.248-273
1244	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1245	!! sur le cycle de l'eau, PhD Thesis, available from:
1246	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1247	!!
1248	!! FLOWCHART : None
1249	!! \n
1250	!_ ================================================================================================================================
1251
1252	SUBROUTINE diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)
1253
1254	!! 0. Variable and parameter declaration
1255
1256	!! 0.1 Input variables
1257	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1258	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
1259	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
1260	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Beta for Interception
1261	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Beta for Transpiration
1262
1263	!! 0.2 Output variables
1264	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
1265
1266	!! 0.3 Modified variables
1267	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! limiting factor for bare soil evaporation
1268	!! when the 11-layer hydrology is used (-)
1269	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! limiting factor for bare soil evaporation
1270	!! when the 11-layer hydrology is used (-)
1271	!! 0.4 Local variables
1272	INTEGER(i_std) :: ji
1273	REAL(r_std), DIMENSION(kjpindex) :: vegtot
1274	!_ ================================================================================================================================
1275
1276	IF (printlev>=3) WRITE (numout,*) 'Entering diffuco_bare'
1277
1278	!! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil
1279
1280	! To use the new version of hydrol_split_soil and ensure the water conservation, we must have throughout sechiba at all time:
1281	! a) evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)soiltile(ji,:)vegtot(ji))
1282	! b) all the terms (vbeta4, evap_bare_lim, evap_bare_lim_ns) =0 if evap_bare_lim(ji) LE min_sechiba
1283	! This must also be kept true in diffuco_comb
1284
1285	DO ji = 1, kjpindex
1286
1287	! The limitation by 1-beta2-beta3 is due to the fact that evaporation under vegetation is possible
1288	!! \latexonly
1289	!! \input{diffucobare3.tex}
1290	!! \endlatexonly
1291
1292	IF ( (evap_bare_lim(ji) .GT. min_sechiba) .AND. &
1293	! in this case we can't have vegtot LE min_sechina, cf hydrol_soil
1294	(un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) .GT. min_sechiba) ) THEN
1295	! eventually, none of the left-hand term is close to zero
1296
1297	vegtot(ji) = SUM(veget_max(ji,:))
1298	IF (evap_bare_lim(ji) < (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)))) THEN
1299	! Standard case
1300	vbeta4(ji) = evap_bare_lim(ji)
1301	ELSE
1302	vbeta4(ji) = un - SUM(vbeta2(ji,:)+vbeta3(ji,:))
1303
1304	! We now have to redefine evap_bare_lim & evap_bare_lim_ns
1305	IF (evap_bare_lim(ji) .GT. min_sechiba) THEN
1306	evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji)
1307	ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles
1308	evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji)
1309	ENDIF
1310
1311	evap_bare_lim(ji) = vbeta4(ji)
1312	! consistent with evap_bare_lim(ji) =
1313	! SUM(evap_bare_lim_ns(ji,:)soiltile(ji,:)vegtot(ji))
1314	! as SUM(soiltile(ji,:)) = 1
1315	END IF
1316
1317	ELSE ! instead of having very small vbeta4, we set everything to zero
1318	vbeta4(ji) = zero
1319	evap_bare_lim(ji) = zero
1320	evap_bare_lim_ns(ji,:) = zero
1321	ENDIF
1322
1323	END DO
1324
1325	IF (printlev>=3) WRITE (numout,*) ' diffuco_bare done '
1326
1327	END SUBROUTINE diffuco_bare
1328
1329
1330	!! ==============================================================================================================================
1331	!! SUBROUTINE : diffuco_trans_co2
1332	!!
1333	!>\BRIEF This subroutine computes carbon assimilation and stomatal
1334	!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987)
1335	!! and Yin X and Struik PC. 2009.
1336	!!
1337	!! DESCRIPTION :\n
1338	!! *** General:\n
1339	!! The equations are different depending on the photosynthesis mode (C3 versus C4).
1340	!! Based on the work of Farquhar, von Caemmerer and Berry who published in 1980 a
1341	!! biochemical model for C3 photosynthetic rates (the FvCB model), Yin and
1342	!! Struik 2009 propose an analytical algorithm that incorporates a gs model and uses gm as a
1343	!! temperature-dependent parameter for calculating assimilation based on the FvCB
1344	!! model for various environmental scenarios. Yin and Struik 2009 also propose a
1345	!! C4-equivalent version of the FvCB model and an associated analytical algorithm
1346	!! for this model. They believe that applying FvCB-type models to both C3 and C4
1347	!! crops is recommended to accurately predict the response of crop photosynthesis
1348	!! to multiple, interactive environmental variables. \n
1349	!!
1350	!! Besides an analytical solution for photosynthesis the scheme includes a
1351	!! modified Arrhenius function for the temperature dependence that accounts for a
1352	!! decrease of Vcmax and Jmax at high temperatures and a temperature dependent
1353	!! J max/Vcmax ratio (Kattge & Knorr 2007). The temperature response of Vcmax and
1354	!! J max was parametrized with values from reanalysed data in literature (Kattge &
1355	!! Knorr 2007), whereas Vcmax and Jmax at a reference temperature of 25Â° C were
1356	!! derived from observed species-specific values in the TRY database (Kattge et
1357	!! al. 2011).
1358	!!
1359	!! In this version here we replace the old trunk assimilation and conductance
1360	!! which were computed over 20 levels of LAI and then integrated at the canopy level.
1361	!! These artificial LAI levels were introduced in order to allow for a saturation of
1362	!! photosynthesis with LAI, but they have no physical meaning. The approach of
1363	! the two stream radiation transfer model was extended for a multi-layer canopy.
1364	!! The user can choose the numbers of layers in the run.def and the available light at
1365	!! each layers decides about the numbers of layers used for the photosynthesis
1366	!! calculation.
1367	!!
1368	!! NOTE: The variable Light_Abs_Tot is the amount of light abosrbed by the canopy,
1369	!! according to the two stream albedo routine. Notice that this is a simple mean
1370	!! between the the absorption profile for light coming from an isotropic source
1371	!! (e.g. clouds, and what light is scattered by aerosols on a sunny day) and the
1372	!! incoming radiation directly from the sun. If the fraction of light from each
1373	!! source is ever provided from the atmospheric model or the forcing file, that
1374	!! will have to be changed in the albedo routines.
1375	!!
1376	!! Also note that I feel sunlight and shaded leaf absorption cannot be calculated
1377	!! by using Collim_Abs_Tot and Isotrop_Abs_Tot from the albedo routines, since
1378	!! direct sunlight will scatter light through the canopy. This light can hit
1379	!! shaded leaves. Therefore, part of Collim_Abs_Tot will be used for sunlight
1380	!! leaves and part for shaded leaves, while all of Isotrop_Abs_Tot will be for
1381	!! shaded leaves.
1382	!!
1383	!! RECENT CHANGE(S):
1384	!!
1385	!! 2018/12 added a nitrogen use efficiency constraint on maximum rate of carboxylation
1386	!!
1387	!! 2018/11 replaced t2m and q2m by temp_air and qair; prognostic variables on first atmospheric model layer
1388	!!
1389	!! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular
1390	!! concentration
1391	!!
1392	!! REFERENCE(S) :
1393	!!
1394	!! Yin X, Struik PC. C3 and C4 photosynthesis models: an overview from the
1395	!! perspective of crop modelling. NJAS-Wageningen Journal of Life Sciences 2009.
1396	!! C.J. Bernacchi, A.R. Portis, H. Nakano, S. von Caemmerer, S.P. Long,
1397	!! Temperature response of mesophyll conductance. Implication for the determination
1398	!! of Rubisco enzyme kinetics and for limitations to photosynthesis in vivo, Plant
1399	!! Physiology 130 (2002) 1992â1998.
1400	!! C.J. Bernacchi, E.L. Singsaas, C. Pimentel, A.R. Portis Jr., S.P. Long,
1401	!! Improved temperature response functions for models of Rubisco-limited
1402	!! photosynthesis, Plant, Cell and Environment 24 (2001) 253â259.
1403	!! B.E. Medlyn, E. Dreyer, D. Ellsworth, M. Forstreuter, P.C. Harley, M.U.F.
1404	!! Kirschbaum, X. Le Roux, P. Montpied, J. Strassemeyer, A. Walcroft, K. Wang, D.
1405	!! Loustau, Temperature response of parameters of a biochemically based model of
1406	!! photosynthesis. II. A review of experimental data, Plant, Cell and Environ- ment
1407	!! 25 (2002) 1167â1179.
1408	!! Kattge, J. and Knorr, W.: Temperature acclimation in a biochemical model of
1409	!! photosynthesis: a reanalysis of data from 36 species, Plant Cell Environ., 30
1410	!! (2007) 1176â1190.
1411	!! Keenan, T., et al., Soil water stress and coupled photosynthesisâconductance
1412	!! models: Bridging the gap between conflicting reports on the relative roles of
1413	!! stomatal, mesophyll conductance and biochemical limitations to photosynthesis.
1414	!! Agric. Forest Meteorol. (2010)
1415	!!
1416	!! FLOWCHART : None
1417	!! \n
1418	!_ ================================================================================================================================
1419
1420	SUBROUTINE diffuco_trans_co2 (kjpindex, swdown, pb, qair, temp_air, &
1421	temp_growth, rau, q_cdrag, humrel, &
1422	assim_param, Ca, lai, &
1423	veget_max, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, &
1424	cimean, gsmean, gpp, cresist, vbeta23, Light_Abs_Tot,&
1425	Light_Tran_Tot, lai_per_level, gs_distribution, gs_diffuco_output, &
1426	gstot_component, gstot_frac, veget, warnings, u_speed, profile_vbeta3, &
1427	profile_rveget, index_assi, leaf_ci, assimtot, assimi_lev, &
1428	hist_id, indexveg, indexlai, index, kjit, cim, info_limitphoto, &
1429	JJ_out)
1430
1431
1432	!
1433	!! 0. Variable and parameter declaration
1434	!
1435
1436	!
1437	!! 0.1 Input variables
1438	!
1439	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1440	REAL(r_std),DIMENSION (:), INTENT (in) :: swdown !! Downwelling short wave flux
1441	!! @tex ($W m^{-2}$) @endtex
1442	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
1443	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Specific humidity at first atmospheric model layer
1444	!! @tex ($kg kg^{-1}$) @endtex
1445	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature at first atmospheric model layer (K)
1446	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
1447	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! air density @tex ($kg m^{-3}$) @endtex
1448	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1449	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (0-1,unitless)
1450	REAL(r_std),DIMENSION (:,:,:), INTENT (in) :: assim_param !! min+max+opt temps (K), vcmax, vjmax for
1451	!! photosynthesis
1452	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1453	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: Ca !! CO2 concentration inside the canopy
1454
1455	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai
1456	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum vegetation fraction of each PFT inside
1457	!! the grid box (0-1, unitless)
1458	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
1459	!! @tex ($kg m^{-2}$) @endte
1460	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
1461	!! @tex ($kg m^{-2}$) @endtex
1462	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will
1463	!! transpire (unitless)
1464	REAL(r_std),DIMENSION (:,:,:),INTENT (in) :: Light_Abs_Tot !! Absorbed radiation per level for photosynthesis
1465	REAL(r_std),DIMENSION (:,:,:),INTENT (in) :: Light_Tran_Tot !! Transmitted radiation per level for photosynthesis
1466	REAL(r_std),DIMENSION(:,:,:),INTENT (in) :: lai_per_level !! This is the LAI per vertical level
1467	!! @tex $(m^{2} m^{-2})$
1468	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT (in) :: veget
1469	REAL(r_std), DIMENSION(:), INTENT (in) :: u_speed !! wind speed at specified canopy level (m/s)
1470
1471	INTEGER(i_std),INTENT(in) :: hist_id !! History_ file identifier (-)
1472	INTEGER(i_std),DIMENSION (kjpindex*nvm),INTENT(in) :: indexveg !! Indices of the points on the 3D map (-)
1473	INTEGER(i_std),DIMENSION (kjpindex*(nlevels_tot+1)),INTENT(in) :: indexlai !! Indices of the points on the 3D map
1474	INTEGER(i_std),DIMENSION (kjpindex),INTENT(in) :: index !! Indices of the points on the map (-)
1475	INTEGER(i_std),INTENT(in) :: kjit !! Time step number (-)
1476
1477	!
1478	!! 0.2 Output variables
1479	!
1480	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration (unitless)
1481	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for Potential Transpiration
1482	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! stomatal resistance of vegetation
1483	!! @tex ($s m^{-1}$) @endtex
1484	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex
1485	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! mean intercellular CO2 concentration
1486	!! @tex ($\mu mol mol^{-1}$) @endtex
1487	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! mean stomatal conductance to CO2 (umol m-2 s-1)
1488	REAL(r_Std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area)
1489	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cresist !! coefficient for resistances (??)
1490
1491	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1492	:: gs_distribution
1493	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1494	:: gs_diffuco_output
1495	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1496	:: gstot_component
1497	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1498	:: gstot_frac
1499	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1500	:: profile_vbeta3
1501	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) &
1502	:: profile_rveget
1503	INTEGER(i_Std),DIMENSION (kjpindex), INTENT (out) :: index_assi !! pixels for a given PFT where photosynthesis takes place
1504	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot), INTENT (out) :: leaf_ci !! intercellular CO2 concentration (ppm)
1505	REAL(r_std),DIMENSION(kjpindex,nvm), INTENT (out) :: assimtot !! total assimilation
1506	REAL(r_std),DIMENSION(kjpindex,nvm,nlevels_tot), INTENT (out) &
1507	:: assimi_lev !! assimilation (at all LAI levels)
1508	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1509
1510	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cim !! Intercellular CO2 over nlevels_tot
1511	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot+1),INTENT(out) &
1512	:: info_limitphoto !! Store information of limitation for photosynthesis
1513
1514	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot),INTENT(out) &
1515	:: JJ_out !! Store JJ
1516	!
1517	!! 0.3 Modified variables
1518	REAL(r_std),DIMENSION(kjpindex,nvm,nwarns), INTENT(inout) &
1519	:: warnings !! A warning counter
1520
1521	!
1522	!! 0.4 Local variables
1523	!
1524
1525	REAL(r_std), DIMENSION(kjpindex,nvm,nlevels_tot+1) :: assimi !! assimilation (at a specific LAI level)
1526	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1527	REAL(r_std), DIMENSION(kjpindex,nvm) :: coupled_lai
1528	!
1529	REAL(r_Std),DIMENSION (kjpindex,nvm,nlevels_tot) :: gs_store
1530	REAL(r_Std),DIMENSION (kjpindex,nvm) :: gs_column_total
1531	REAL(r_Std),DIMENSION (kjpindex,nvm) :: gstot_store
1532
1533
1534	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_gs_store
1535	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_cimean
1536	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_assimtot
1537
1538	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_Rdtot
1539	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_gstot
1540	REAL(r_std),DIMENSION (kjpindex,nlevels_tot) :: profile_gstop
1541
1542	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_gsmean
1543
1544	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_gpp
1545
1546
1547	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: profile_rstruct
1548	REAL(r_std),DIMENSION (nlevels_tot) :: profile_speed
1549	REAL(r_std),DIMENSION (nlevels_tot) :: profile_cresist
1550
1551	!
1552	!
1553	REAL(r_std),DIMENSION(kjpindex,nvm,nlevels_tot) :: Abs_light !! fraction of absorbed light within each layer (just to check)
1554	REAL(r_std),DIMENSION(kjpindex,nvm) :: Abs_light_can !! fraction of absorbed light within the canopy (just to check)
1555	REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax !! maximum rate of carboxylation Notice that,
1556	!! as it is currently written, this is only
1557	!! a function of nue and leafN_top
1558	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1559	REAL(r_std) :: N_profile !! N-content reduction factor for a specific canopy depth (unitless)
1560	REAL(r_std) :: N_normalisation !! Normalisation factor for the N-profile,
1561	!! to get the correct total N content (unitless)
1562	REAL(r_std),DIMENSION (kjpindex) :: leafN_top !! Leaf N content - top of canopy (gN m-2[leaf])
1563	INTEGER(i_std) :: ji, jv, ilev, jl, limit_photo !! indices (unitless)
1564	REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest !! intercellular CO2 concentration at the lowest
1565	!! LAI level
1566	!! @tex ($\mu mol mol^{-1}$) @endtex
1567	REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap !! relative water quantity in the water
1568	!! interception reservoir (0-1,unitless)
1569	REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex
1570	! Assimilation
1571	LOGICAL, DIMENSION(kjpindex) :: assimilate !! where assimilation is to be calculated
1572	!! (unitless)
1573	INTEGER(i_std) :: nia,inia,nina !! counter/indices (unitless)
1574	INTEGER(i_std) :: inina,iainia !! counter/indices (unitless)
1575	INTEGER(i_std), DIMENSION(kjpindex) :: index_non_assi !! indices (unitless)
1576	REAL(r_std), DIMENSION(kjpindex, nlevels_tot+1) :: vc2 !! rate of carboxylation (at a specific LAI level)
1577	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1578
1579	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: vj2 !! rate of Rubisco regeneration (at a specific LAI
1580	!! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex
1581	REAL(r_std), DIMENSION(kjpindex) :: gstop !! stomatal conductance to H2O at topmost level
1582	!! @tex ($m s^{-1}$) @endtex
1583
1584	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: gs !! stomatal conductance to CO2
1585	!! @tex ($\mol m^{-2} s^{-1}$) @endtex
1586	REAL(r_std), DIMENSION(kjpindex) :: gamma_star !! CO2 compensation point (ppm)
1587	!! @tex ($\mu mol mol^{-1}$) @endtex
1588	REAL(r_std), DIMENSION(kjpindex) :: air_relhum !! air relative humidity at 2m
1589	!! @tex ($kg kg^{-1}$) @endtex
1590	REAL(r_std), DIMENSION(kjpindex) :: VPD !! Vapor Pressure Deficit (kPa)
1591	REAL(r_std), DIMENSION(kjpindex) :: water_lim !! water limitation factor (0-1,unitless)
1592
1593	REAL(r_std), DIMENSION(kjpindex) :: gstot !! total stomatal conductance to H2O
1594	!! Final unit is
1595	!! @tex ($m s^{-1}$) @endtex
1596	REAL(r_std), DIMENSION(kjpindex) :: Rdtot !! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1)
1597	REAL(r_std), DIMENSION(kjpindex) :: gs_top !! leaf stomatal conductance to H2O across all top levels
1598	!! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex
1599	REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??)
1600	!! @tex ($g g^{-1}$) @endtex
1601	!! levels (0-1,unitless)
1602	REAL(r_std), DIMENSION(kjpindex) :: T_Vcmax !! Temperature dependance of Vcmax (unitless)
1603	REAL(r_std), DIMENSION(kjpindex) :: S_Vcmax_acclim_temp !! Entropy term for Vcmax
1604	!! accounting for acclimation to temperature (J K-1 mol-1)
1605	REAL(r_std), DIMENSION(kjpindex) :: T_Jmax !! Temperature dependance of Jmax
1606	REAL(r_std), DIMENSION(kjpindex) :: S_Jmax_acclim_temp !! Entropy term for Jmax
1607	!! accounting for acclimation toxs temperature (J K-1 mol-1)
1608	REAL(r_std), DIMENSION(kjpindex) :: T_gm !! Temperature dependance of gmw
1609	REAL(r_std), DIMENSION(kjpindex) :: T_Rd !! Temperature dependance of Rd (unitless)
1610	REAL(r_std), DIMENSION(kjpindex) :: T_Kmc !! Temperature dependance of KmC (unitless)
1611	REAL(r_std), DIMENSION(kjpindex) :: T_KmO !! Temperature dependance of KmO (unitless)
1612	REAL(r_std), DIMENSION(kjpindex) :: T_Sco !! Temperature dependance of Sco
1613	REAL(r_std), DIMENSION(kjpindex) :: T_gamma_star !! Temperature dependance of gamma_star (unitless)
1614	REAL(r_std), DIMENSION(kjpindex) :: vc !! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 mâ2 sâ1)
1615	REAL(r_std), DIMENSION(kjpindex) :: vj !! Maximum rate of e- transport under saturated light (mumol CO2 mâ2 sâ1)
1616	REAL(r_std), DIMENSION(kjpindex) :: gm !! Mesophyll diffusion conductance (molCO2 mÃ¢2 sÃ¢1 barÃ¢1)
1617	REAL(r_std), DIMENSION(kjpindex) :: g0var
1618	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: Rd !! Day respiration (respiratory CO2 release other than by photorespiration)(mumol CO2 mÃ¢2 sÃ¢1)
1619	REAL(r_std), DIMENSION(kjpindex) :: Kmc !! MichaelisâMenten constant of Rubisco for CO2 (mubar)
1620	REAL(r_std), DIMENSION(kjpindex) :: KmO !! MichaelisâMenten constant of Rubisco for O2 (mubar)
1621	REAL(r_std), DIMENSION(kjpindex) :: Sco !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1)
1622	REAL(r_std), DIMENSION(kjpindex) :: gb_co2 !! Boundary-layer conductance (molCO2 mÃ¢2 sÃ¢1 barÃ¢1)
1623	REAL(r_std), DIMENSION(kjpindex) :: gb_h2o !! Boundary-layer conductance (molH2O mÃ¢2 sÃ¢1 barÃ¢1)
1624	REAL(r_std), DIMENSION(kjpindex) :: fvpd !! Factor for describing the effect of leaf-to-air vapour difference on gs (-)
1625	REAL(r_std), DIMENSION(kjpindex) :: low_gamma_star !! Half of the reciprocal of Sc/o (bar bar-1)
1626	REAL(r_std) :: fcyc !! Fraction of electrons at PSI that follow cyclic transport around PSI (-)
1627	REAL(r_std) :: z !! A lumped parameter (see Yin et al. 2009) ( mol mol-1)
1628	REAL(r_std) :: Rm !! Day respiration in the mesophyll (umol CO2 mâ2 sâ1)
1629	REAL(r_std) :: Cs_star !! Cs -based CO2 compensation point in the absence of Rd (ubar)
1630	REAL(r_std), DIMENSION(kjpindex) :: Iabs !! Photon flux density absorbed by leaf photosynthetic pigments (umol photon mâ2 sâ1)
1631	REAL(r_std), DIMENSION(kjpindex) :: Jmax !! Maximum value of J under saturated light (umol eâ mâ2 sâ1)
1632	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: JJ !! Rate of eâ transport (umol eâ mâ2 sâ1)
1633	REAL(r_std) :: J2 !! Rate of all eâ transport through PSII (umol eâ mâ2 sâ1)
1634	REAL(r_std) :: VpJ2 !! eâ transport-limited PEP carboxylation rate (umol CO2 mâ2 sâ1)
1635	REAL(r_std) :: A_1, A_3 !! Lowest First and third roots of the analytical solution for a general
1636	!! cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
1637	REAL(r_std) :: A_1_tmp, A_3_tmp !! Temporary First and third roots of the analytical solution
1638	!! for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
1639	REAL(r_std) :: Obs !! Bundle-sheath oxygen partial pressure (ubar)
1640	REAL(r_std), DIMENSION(kjpindex,nlevels_tot+1) :: Cc_loc !! Chloroplast CO2 partial pressure (ubar)
1641	!! Notice we had _loc to the name to avoid
1642	!! conflict with a variable already declared
1643	!! in constantes_var
1644	REAL(r_std) :: ci_star !! Ci -based CO2 compensation point in the absence of Rd (ubar)
1645	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))
1646	REAL(r_std) :: QQ,UU,PSI,x1,x2,x3 !! Variables used for solving the cubic equation (see Yin et al. (2009))
1647	rEAL(r_std) :: lai_min !! mininum lai to calculate photosynthesis
1648	REAL(r_std), DIMENSION(kjpindex) :: laisum !! when calculating cim over nlevels_tot
1649
1650	LOGICAL,DIMENSION(kjpindex,nlevels_tot) :: top_level !! A flag to see if a given level is considered
1651	!! a "top" level or not. Must be done this way
1652	!! because there could be several "top" levels.
1653	REAL(r_std) :: totlai_top !! The sum of the LAI in all the levels that we
1654	!! declare as the "top".
1655	INTEGER(i_std) :: lev_count !! How many levels we have already taken for
1656	!! the "top"
1657
1658	! @defgroup Photosynthesis Photosynthesis
1659	! @{
1660	! 1. Preliminary calculations\n
1661	REAL(r_std),DIMENSION (kjpindex,nvm,nlevels_tot) :: light_tran_to_level !! Cumulative amount of light transmitted to a given level
1662	LOGICAL,DIMENSION (kjpindex) :: lskip !! A flag to indicate we've hit an unusual
1663	!! case were we may have numerical instability.
1664	INTEGER :: ic, jc, kc
1665
1666
1667	INTEGER :: ipts
1668	LOGICAL :: llight
1669	LOGICAL :: llai
1670	INTEGER,SAVE :: istep=0
1671	!$OMP THREADPRIVATE(istep)
1672
1673	!_ ================================================================================================================================
1674
1675	leaf_ci(:,:,:) = zero
1676
1677	IF (printlev>=4) WRITE (numout,*) ' diffuco_trans_co2'
1678	IF (printlev_loc>=5) WRITE(numout,*) 'diffuco_trans_co2: nlevels_tot,lai_per_level',nlevels_tot,lai_per_level(:,:,:)
1679
1680	gs_store = 0.0d0
1681	gs_column_total = 0.0d0
1682	gs_distribution = 0.0d0
1683	gs_diffuco_output = 0.0d0
1684
1685	gstot_frac = 0.0d0
1686	gstot_component = 0.d0
1687	gstot_store = 0.0d0
1688
1689	profile_gs_store = 0.0d0
1690	profile_cimean = 0.0d0
1691	profile_assimtot = 0.0d0
1692
1693	profile_Rdtot = 0.0d0
1694	profile_gstot = 0.0d0
1695	profile_gstop = 0.0d0
1696
1697	profile_gsmean = 0.0d0
1698	profile_gpp = 0.0d0
1699	profile_rveget = 2.5d20
1700	profile_rstruct = 0.0d0
1701	profile_speed = 0.0d0
1702	profile_cresist = 0.0d0
1703	profile_vbeta3 = 0.0d0
1704
1705	JJ_out = 0.0d0
1706	! Debug
1707	IF(printlev_loc>=4)THEN
1708	WRITE(numout,*) 'Absorbed light in diffuco. ',test_grid,test_pft
1709	DO ilev=nlevels_tot,1,-1
1710	WRITE(numout,'(I5,10F20.10)') &
1711	ilev, Light_Abs_Tot(test_grid,test_pft,ilev),&
1712	Light_Tran_Tot(test_grid,test_pft,ilev),&
1713	lai_per_level(test_grid,test_pft,ilev)
1714	ENDDO
1715	ENDIF
1716	!-
1717
1718	gs_column_total = 0.0d0
1719	gs_distribution = 0.0d0
1720	gs_diffuco_output = 0.0d0
1721	gstot_frac = 0.0d0
1722	gstot_component = 0.d0
1723	gstot_store = 0.0d0
1724	gs_store(:,:,:)=zero
1725	!
1726	! @addtogroup Photosynthesis
1727	! @{
1728	! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance).\n
1729	!! \latexonly
1730	!! \input{diffuco_trans_co2_1.3.tex}
1731	!! \endlatexonly
1732	! @}
1733	!
1734	CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
1735	air_relhum(:) = &
1736	( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) / &
1737	( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) )
1738
1739	! WRITE(numout, *) '050315 air_relhum is: ', air_relhum
1740
1741	VPD(:) = ( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) ) &
1742	- ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) )
1743	! VPD is needed in kPa
1744	! We limit the impact of VPD in the range of [0:6] kPa
1745	! This seems to be the typical range of values
1746	! that the vapor pressure difference can take, although
1747	! there is a paper by Mott and Parkhurst, Plant, Cell & Environment,
1748	! Vol. 14, pp. 509--519 (1991) which gives a value of 15 kPa.
1749	!!$ VPD(:) = MAX(zero,MIN(VPD(:)/10.,6.))
1750	VPD(:) = VPD(:)/10.
1751
1752	!
1753	! 2. beta coefficient for vegetation transpiration
1754	!
1755	rstruct(:,1) = rstruct_const(1)
1756	!+++CHECK+++
1757	! Setting rveget makes the model crash when running in debug mode becaus
1758	! the default xios settings cannot write undefined values.
1759	!!$ rveget(:,:) = undef_sechiba
1760	rveget(:,:) = un
1761	!+++++++++++
1762	vbeta3(:,:) = zero
1763	vbeta3pot(:,:) = zero
1764	gsmean(:,:) = zero
1765	gpp(:,:) = zero
1766	light_tran_to_level(:,:,:)=un
1767	assimi_lev(:,:,:) = zero
1768	cim(:,:) = zero
1769	cimean(:,1) = Ca(:)
1770	info_limitphoto(:,:,:) = zero
1771	!
1772	! 2.1 Initializations
1773	! Loop over vegetation types
1774	!
1775
1776	DO jv = 2,nvm
1777	gamma_star(:) = zero
1778	Kmo(:) = zero
1779	Kmc(:) = zero
1780	gm(:) = zero
1781	g0var(:) = zero
1782
1783	Cc_loc(:,:) = zero
1784	vc2(:,:) = zero
1785	vj2(:,:) = zero
1786	JJ(:,:) = zero
1787	gs(:,:) = zero
1788	assimi(:,:,:) = zero
1789	Rd(:,:) = zero
1790
1791	!
1792	! beta coefficient for vegetation transpiration
1793	!
1794	rstruct(:,jv) = rstruct_const(jv)
1795	cimean(:,jv) = Ca(:)
1796
1797	! ---TEMP---
1798	! calculation of vcmax before N is included
1799	! IF (downregulation_co2) THEN
1800	! vcmax(:,jv) = assim_param(:,jv,ivcmax)*&
1801	! (un-downregulation_co2_coeff(jv)*&
1802	! log(Ca(:)/downregulation_co2_baselevel))
1803	! ELSE
1804	! vcmax(:,jv) = assim_param(:,jv,ivcmax)
1805	! ENDIF
1806	!-----------
1807
1808	WHERE (lai(:,jv) .GT. min_sechiba)
1809	!The normalisation factor is now computed later, based on
1810	!the profile of the transmitted light.
1811	!leafN_top(:) = assim_param(:,jv,ileafN) * ext_coeff_N(jv) / &
1812	! ( 1. -exp(-ext_coeff_N(jv) * lai(:,jv)) )
1813	leafN_top(:) = assim_param(:,jv,ileafN)
1814
1815	ELSEWHERE
1816
1817	leafN_top(:) = zero
1818
1819	ENDWHERE
1820
1821	! Very import line of code. Here vcmas is CALCULATED from nue and ileafN
1822	! which are both stored in assim_param. Note that the calculation of
1823	! vcmax in stomate_vmax is no longer used in the model.
1824	! Note here leaf nitrogen per ground area is used (by integrating all leaves),
1825	! thus vcmax is calculated based on ground area as well. The vcmax will be
1826	! converted to be based per leaf area per LAI layer later on.
1827	vcmax(:,jv) = assim_param(:,jv,inue)*leafN_top(:)
1828
1829	!WRITE(numout,*) '210515 assim_param(grid,pft,npco2)','grid:'grid,'pft:',pft,'npco2:',npco2
1830	!WRITE(numout,*) assim_param(:,:,:)
1831	!! mask that contains points where there is photosynthesis
1832	!! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points.
1833	!! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not
1834	!! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation,
1835	!! temperature and relative humidity).
1836	!! For the points where assimilation is not calculated, variables are initialized to specific values.
1837	!! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
1838	!! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes
1839	!! the nina points with no assimilation.
1840	!
1841	! set minimum lai value to test if photosynthesis should be calculated
1842	! 0.0001 is abitrary
1843	! could be externalised
1844	lai_min = 0.0001
1845	!
1846	nia=0
1847	nina=0
1848	index_assi=0
1849	index_non_assi=0
1850	llai=.FALSE.
1851	!
1852	DO ji=1,kjpindex
1853
1854	IF(veget_max(ji,jv) == zero)THEN
1855	! this vegetation type is not present, so no reason to do the
1856	! calculation. However, WHERE loops later on depend on the
1857	! value of assimilate, so we set it here.
1858	assimilate(ji) = .FALSE.
1859	CYCLE
1860	ENDIF
1861
1862	!
1863	! This is tricky. We need a certain amount of LAI to have photosynthesis.
1864	! However, we divide our canopy up into certain levels, so even if we
1865	! have a total amount that is sufficient, there might not be enough
1866	! in any particular layer. So we check each layer to see if there
1867	! is one that has enough to give a non-zero photosynthesis. In
1868	! the worst-case scenario (small LAI in every other level), this becomes
1869	! our top level for the gs calculation. If gs_top_level is zero below,
1870	! we have a divide by zero issue.
1871	llai=.FALSE.
1872	DO ilev=1,nlevels_tot
1873	IF(lai_per_level(ji,jv,ilev) .GT. lai_min) THEN
1874	llai=.TRUE.
1875	ELSE
1876	ENDIF
1877	ENDDO
1878
1879
1880	IF ( llai .AND. &
1881	( veget_max(ji,jv) .GT. min_sechiba ) ) THEN
1882
1883	! Here we need to do something different than in the previous model.
1884	! We explicitly need to check for absorbed light. This occurs
1885	! because our albedo in the two stream model is a function
1886	! of the solar angle. If the sun has set, we have no
1887	! light, and therefore we cannot have photosynthesis.
1888	llight=.FALSE.
1889	DO ilev=1,nlevels_tot
1890	IF(Light_Abs_Tot(ji,jv,ilev) .GT. min_sechiba) llight=.TRUE.
1891	ENDDO
1892	IF ( ( swdown(ji) .GT. min_sechiba ) .AND. &
1893	( humrel(ji,jv) .GT. min_sechiba) .AND. &
1894	( temp_growth(ji) .GT. tphoto_min(jv) ) .AND. &
1895	llight .AND. &
1896	( temp_growth(ji) .LT. tphoto_max(jv) )) THEN
1897	!
1898	assimilate(ji) = .TRUE.
1899	nia=nia+1
1900	index_assi(nia)=ji
1901	!
1902	ELSE
1903	!
1904	assimilate(ji) = .FALSE.
1905	nina=nina+1
1906	index_non_assi(nina)=ji
1907	!
1908	ENDIF
1909	ELSE
1910	!
1911	assimilate(ji) = .FALSE.
1912	nina=nina+1
1913	index_non_assi(nina)=ji
1914	!
1915	ENDIF
1916	!
1917	ENDDO
1918	!
1919
1920	IF(printlev_loc>=4 .AND. jv == test_pft)THEN
1921	WRITE(numout,*) 'jv: ',jv
1922	WRITE(numout,*) 'test_grid: ',test_grid
1923	WRITE(numout,*) 'llai, veget_max: ',llai,veget_max(:,jv)
1924	WRITE(numout,*) 'temp_growth: ',temp_growth(test_grid)
1925	WRITE(numout,*) 'tphoto_min: ',tphoto_min(jv)
1926	WRITE(numout,*) 'tphoto_max: ',tphoto_max(jv)
1927	WRITE(numout,*) 'swdown: ',swdown(test_grid)
1928	WRITE(numout,*) 'index_non_assi(:): ',index_non_assi(test_grid)
1929	WRITE(numout,*) 'index_assi(:): ',index_assi(test_grid)
1930	WRITE(numout,*) 'nia: ',nia
1931	WRITE(numout,*) 'nina',nina
1932	ENDIF
1933
1934	top_level(:,:)=.FALSE.
1935	gstot(:) = zero
1936	gstop(:) = zero
1937	assimtot(:,jv) = zero
1938	Rdtot(:)=zero
1939	gs_top(:) = zero
1940	lskip(:)=.FALSE.
1941	!
1942	zqsvegrap(:) = zero
1943	WHERE (qsintmax(:,jv) .GT. min_sechiba)
1944	!! relative water quantity in the water interception reservoir
1945	zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
1946	ENDWHERE
1947
1948	!! Calculate the water limitation factor that restricts Vcmax when hydrol_arch is turned off.
1949	! If hydrol_architecture is used, VCmax is not controled but instead there is a control on
1950	! stomatal conductance
1951	IF(ok_hydrol_arch) THEN
1952	WHERE ( assimilate(:) )
1953	water_lim(:) = un
1954	ELSEWHERE
1955	water_lim(:) = zero
1956	ENDWHERE
1957	ELSE
1958	WHERE ( assimilate(:) )
1959	water_lim(:) = humrel(:,jv)
1960	ELSEWHERE
1961	water_lim(:) = zero
1962	ENDWHERE
1963	ENDIF ! (ok_hydrol_arch)
1964
1965	! give a default value of ci for all pixel that do not assimilate
1966	DO ilev=1,nlevels_tot
1967	DO inina=1,nina
1968	leaf_ci(index_non_assi(inina),jv,ilev) = Ca(index_non_assi(inina))
1969	ENDDO
1970	ENDDO
1971
1972	!
1973	!
1974	! Here is the calculation of assimilation and stomatal conductance
1975	! based on the work of Farquahr, von Caemmerer and Berry (FvCB model)
1976	! as described in Yin et al. 2009
1977	! Yin et al. developed a extended version of the FvCB model for C4 plants
1978	! and proposed an analytical solution for both photosynthesis pathways (C3 and C4)
1979	! Photosynthetic parameters used are those reported in Yin et al.
1980	! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use
1981	! Medlyn et al. (2002) and Kattge & Knorr (2007)
1982	! Because these 2 references do not consider mesophyll conductance, we neglect this term
1983	! in the formulations developed by Yin et al.
1984	! Consequently, gm (the mesophyll conductance) tends to the infinite
1985	! This is of importance because as stated by Kattge & Knorr and Medlyn et al.,
1986	! values of Vcmax and Jmax derived with different model parametrizations are not
1987	! directly comparable and the published values of Vcmax and Jmax had to be standardized
1988	! to one consistent formulation and parametrization
1989
1990	! See eq. 6 of Yin et al. (2009)
1991	! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
1992	T_KmC(:) = Arrhenius(kjpindex,temp_air,298.,E_KmC(jv))
1993	T_KmO(:) = Arrhenius(kjpindex,temp_air,298.,E_KmO(jv))
1994	T_Sco(:) = Arrhenius(kjpindex,temp_air,298.,E_Sco(jv))
1995	T_gamma_star(:) = Arrhenius(kjpindex,temp_air,298.,E_gamma_star(jv))
1996
1997
1998	! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
1999	T_Rd(:) = Arrhenius(kjpindex,temp_air,298.,E_Rd(jv))
2000
2001
2002	! For C3 plants, we assume that the Entropy term for Vcmax and Jmax
2003	! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10
2004	! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function
2005	! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002)
2006	! and Kattge & Knorr (2007).
2007	! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function
2008	! Concerning this apparent unconsistency, have a look to the section 'Limitation of
2009	! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation
2010
2011	! Growth temperature tested by Kattge & Knorr range from 11 to 35Â°C
2012	! So, we limit the relationship between these lower and upper limits
2013	S_Jmax_acclim_temp(:) = aSJ(jv) + bSJ(jv) * MAX(11., MIN(temp_growth(:),35.))
2014	T_Jmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Jmax(jv),D_Jmax(jv),S_Jmax_acclim_temp)
2015
2016	S_Vcmax_acclim_temp(:) = aSV(jv) + bSV(jv) * MAX(11., MIN(temp_growth(:),35.))
2017	T_Vcmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
2018
2019	vc(:) = vcmax(:,jv) * T_Vcmax(:)
2020
2021
2022	! As shown by Kattge & Knorr (2007), we make use
2023	! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants
2024	! rJV is written as a function of the growth temperature
2025	! rJV = arJV + brJV * T_month
2026	! See eq. 10 of Kattge & Knorr (2007)
2027	! and Table 3 for Values of arJV anf brJV
2028	! Growth temperature is monthly temperature (expressed in Â°C) - See first paragraph of
2029	! section Methods/Data of Kattge & Knorr
2030	vj(:) = ( arJV(jv) + brJV(jv) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:,jv) * T_Jmax(:)
2031
2032	T_gm(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_gm(jv),D_gm(jv),S_gm(jv))
2033	gm(:) = gm25(jv) * T_gm(:) * MAX(1-stress_gm(jv), water_lim(:))
2034
2035	g0var(:) = g0(jv)* MAX(1-stress_gs(jv), water_lim(:))
2036	! @endcodeinc
2037	!
2038	KmC(:)=KmC25(jv)*T_KmC(:)
2039	KmO(:)=KmO25(jv)*T_KmO(:)
2040	Sco(:)=Sco25(jv)*T_sco(:)
2041	gamma_star(:) = gamma_star25(jv)*T_gamma_star(:)
2042	! low_gamma_star is defined by Yin et al. (2009)
2043	! as the half of the reciprocal of Sco - See Table 2
2044	!!$ We derive its value from Equation 7 of Yin et al. (2009)
2045	!!$ assuming a constant O2 concentration (Oi)
2046	!!$ low_gamma_star(:) = gamma_star(:)/Oi
2047	low_gamma_star(:) = 0.5 / Sco(:)
2048
2049	! VPD expressed in kPa
2050	! Note : MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:))))
2051	! is always between 0-1 not including 0 and 1
2052	fvpd(:) = 1. / ( 1. / MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) - 1. ) &
2053	* MAX(1-stress_gs(jv), water_lim(:))
2054
2055	! leaf boundary layer conductance
2056	! conversion from a conductance in (m s-1) to (mol H2O m-2 s-1)
2057	! from Pearcy et al. (1991, see below)
2058	gb_h2o(:) = gb_ref * 44.6 * (tp_00/temp_air(:)) * (pb(:)/pb_std)
2059
2060	! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1)
2061	gb_co2(:) = gb_h2o(:) / ratio_H2O_to_CO2
2062
2063	!
2064	! 2.4 Loop over LAI levels to estimate assimilation and conductance
2065	!
2066	JJ(:,:)=zero
2067	vc2(:,:)=zero
2068	vj2(:,:)=zero
2069	Cc_loc(:,:)=zero
2070	gs(:,:)=zero
2071	Rd(:,:)=zero
2072
2073	IF (nia .GT. 0) THEN
2074
2075	assim_loop:DO inia=1,nia
2076
2077	!
2078	iainia=index_assi(inia)
2079
2080	!
2081	! Find the top layer that still contains lai. We have also
2082	! added a check in here to make sure this layer has absorbed
2083	! light. If there is no absorbed light, we will not calculate
2084	! stomatal conductance for the top level, which means gstop
2085	! will be zero later and we'll crash when we divide by it.
2086	! Include a warning so we know how often this happens.
2087	! The old system of chosing just one top level caused some bad results
2088	! in the transpir, evidently due to the calculation of vbeta3 being
2089	! too low due to too low of LAI in the top level. Therefore, we are
2090	! now going to take several top levels in an effort to make sure
2091	! we have enough LAI. We aim for an arbitrary amount of 0.1 or
2092	! the top three levels. This has not been rigoursly tested, but it
2093	! is theoretically justified by the fact that if the LAI of the top
2094	! level is very low, other leaves will be getting a lot of sun and
2095	! well-exposed to the atmosphere, and therefore contribute equal
2096	! to the resistence of the vegetation.
2097	totlai_top=zero
2098	lev_count=0
2099	!
2100	DO ilev=nlevels_tot, 1, -1
2101	! lai_min is 0.0001. If a layer with so little lai occurs at the
2102	! bottom of the canopy, the absorbed light can be below 10e-8.
2103	! This implies that we should test wether Light_Abs_Tot exceeds
2104	! zero rather than 10-8 (min_sechiba). It is difficult to set a
2105	! strict threshold above which Light_Abs_Tot should be if
2106	! the lai_per_level is above lai_min. Setting the threshold to
2107	! zero is save but not very useful as it will only track the
2108	! worst possible errors.
2109	IF(lai_per_level(iainia,jv,ilev) .GT. lai_min)THEN
2110	IF(Light_Abs_Tot(iainia,jv,ilev) .GE. zero)THEN
2111	IF(totlai_top .LT. lai_top(jv) .AND. lev_count .LT. nlev_top)THEN
2112	lev_count=lev_count+1
2113	top_level(iainia,ilev)=.TRUE.
2114	totlai_top=totlai_top+lai_per_level(iainia,jv,ilev)
2115	ELSE
2116	! We have enough levels now.
2117	!WRITE(numout,*) 'Top levels ',totlai_top,lev_count
2118	EXIT
2119	ENDIF
2120	ELSE
2121	WRITE(numout,*) 'ERROR: Found a level with LAI, ' // &
2122	'but no absorbed light.'
2123	WRITE(numout,*) 'ERROR: iainia,jv,ilev, ',iainia,jv,ilev
2124	WRITE(numout,*) 'ERROR: Light_Abs_Tot(iainia,jv,ilev),' //&
2125	'lai_per_level(iainia,jv,ilev),',iainia,jv, &
2126	ilev,Light_Abs_Tot(iainia,jv,ilev), &
2127	lai_per_level(iainia,jv,ilev)
2128	CALL ipslerr_p(3,'diffuco_trans_co2', &
2129	'Found a level with LAI but no absorbed light','','')
2130	ENDIF
2131	ENDIF
2132	ENDDO
2133	IF(totlai_top .LT. min_sechiba)THEN
2134	! There seem to be some fringe cases where this happens, usually
2135	! at sunrise or sunset with low LAI levels. Instead of stopping
2136	! here, I'm going to issue a warning. If this case happens,
2137	! we produce no GPP, but we're not going to crash. We need to
2138	! be careful in one of the loops below, since a stomatal
2139	! conductance of zero will cause a crash.
2140	IF(err_act.GT.1)THEN
2141	WRITE(numout,*) 'WARNING: Did not find a level with LAI in diffuco!'
2142	WRITE(numout,*) 'WARNING: iainia,ivm: ',iainia,jv
2143	ENDIF
2144	assimtot(iainia,jv) = zero
2145	Rdtot(iainia) = zero
2146	gstot(iainia) = zero
2147	assimi_lev = zero
2148	lskip(iainia)=.TRUE.
2149	CYCLE assim_loop
2150	ENDIF
2151
2152	! We need to know how much is transmitted to this level
2153	! from the top of the canopy, since this absolute number
2154	! is what is actually used in the photosyntheis. Right now
2155	! Light_Tran_Tot and Light_Abs_Tot are the
2156	! relative quantities transmitted through and absorbed
2157	! by the layer.
2158	DO ilev = nlevels_tot-1,1,-1
2159
2160	light_tran_to_level(iainia,jv,ilev) = &
2161	light_tran_to_level(iainia,jv,ilev+1) * &
2162	Light_Tran_Tot(iainia,jv,ilev+1)
2163
2164	ENDDO
2165
2166	!We sum N_profile over the levels where it will be used in the next loop below.
2167	!The total will give the normalisation factor.
2168	N_normalisation=0
2169	DO ilev = 1, nlevels_tot
2170	IF (Light_Abs_Tot(iainia,jv,ilev) .GT. min_sechiba) THEN
2171	IF(lai_per_level(iainia,jv,ilev) .LE. min_sechiba)THEN
2172	CYCLE
2173	ENDIF
2174	N_profile = ( un - .7_r_std * &
2175	( un - light_tran_to_level(iainia,jv,ilev) ) )
2176	N_normalisation = N_normalisation+N_profile*lai_per_level(iainia,jv,ilev)
2177	ENDIF
2178	ENDDO
2179
2180	! Debug
2181	IF(printlev_loc>=4 .AND. jv == test_pft .AND. iainia == test_grid)THEN
2182	WRITE(numout,*) 'Absolute transmitted light to each level',&
2183	jv,iainia
2184	DO ilev = nlevels_tot,1,-1
2185	WRITE(numout,*) ilev,light_tran_to_level(iainia,jv,ilev)
2186
2187	ENDDO
2188	ENDIF
2189	!-
2190
2191	DO ilev = 1, nlevels_tot
2192
2193	IF (Light_Abs_Tot(iainia,jv,ilev) .GT. min_sechiba) THEN
2194
2195	IF(lai_per_level(iainia,jv,ilev) .LE. min_sechiba)THEN
2196
2197	! This is not necessarily a problem. LAI is updated in
2198	! slowproc, while absorbed light is updated in condveg,
2199	! which is called after diffuco. So we might be one
2200	! timestep off. It generally should be at night,
2201	! though, so there is no absorbed light. If this
2202	! persists, that is a real problem.
2203	IF(err_act.GT.1)THEN
2204	WRITE(numout,*) 'WARNING: We have absorbed light but no LAI ' // &
2205	'in this level. Skipping.'
2206	WRITE(numout,'(A,3I8)') 'WARNING: iainia,jv,ilev: ',iainia,jv,ilev
2207	WRITE(numout,'(A,2F20.10)') 'WARNING: Light_Abs_Tot(iainia,jv,ilev), ' // &
2208	'lai_per_level(iainia,jv,ilev): ',&
2209	Light_Abs_Tot(iainia,jv,ilev), &
2210	lai_per_level(iainia,jv,ilev)
2211	ENDIF
2212	CYCLE
2213	ENDIF
2214
2215	! nitrogen is scaled according to the transmitted light
2216	! before the caclulation was: N_profile =
2217	! exp( -ext_coeff_N(jv)*laitab(jl) ). Where laitab is the
2218	! lai contained in a fixed lai layer. ext_coeff_N is
2219	! prescribed so this was a prescribed profile. We replaced
2220	! it by the transmitted light calculated by the 2stream
2221	! model, maybe ( un - Light_Tran_Tot(iainia,jv,ilev) )
2222	! could be replaced by the directly calculated absorbed
2223	! light Light_Abs_Tot(iainia,jv,ilev).
2224	! Notice that light in the old code was the cumulative
2225	! light transmitted to the layer, while Light_Tran_Tot
2226	! is the relative light transmitted through the layer.
2227	! Note the trick to convert vcmax from per ground area to
2228	! per leaf area is made here and below. N_profile calculated here
2229	! and used below, is implicitly divided by the leaf area index
2230	! (leaf area per ground area). c.f. how N_normalisation is calculated
2231	! above.
2232	N_profile = ( un - .7_r_std * &
2233	( un - light_tran_to_level(iainia,jv,ilev) ) )/N_normalisation
2234
2235	! the vcmax and photosynthesis is only restricted by
2236	! stomatal closure when transpiration is higher than the
2237	! water supply to the leaves as calculated in hydrol_arch.
2238	! This more consistent with the substantial consensus in
2239	! plant physiology that the main restriction of
2240	! photosynthesis during drought occurs through diffusion
2241	! limitation rather than photochemical or biochemical failure
2242	! (Flexas et al. 2006). Effects of drought stress on
2243	! respiration are less well known: increases, decreases and
2244	! no change have been found. For now there is no effect of
2245	! drought stress on dark respiration when hydrol_arch is
2246	! used.
2247	vc2(iainia,ilev) = vc(iainia) * N_profile * MAX(1-stress_vcmax(jv), water_lim(iainia))
2248	vj2(iainia,ilev) = vj(iainia) * N_profile * MAX(1-stress_vcmax(jv), water_lim(iainia))
2249
2250	! see Comment in legend of Fig. 6 of Yin et al. (2009)
2251	! Rd25 is assumed to equal 0.01 Vcmax25
2252	Rd(iainia,ilev) = vcmax(iainia,jv) * N_profile * 0.01 * &
2253	T_Rd(iainia) * MAX(1-stress_vcmax(jv), water_lim(iainia))
2254	! As Light_Abs_Tot is the absorbed radiation per layer,
2255	! expressed as the fraction of overall light hitting the
2256	! canopy, there is no need for an additional multiplication
2257	! by light_tran_to_level, which is is the cumulative
2258	! fraction of transmitted light (otherwise we are counting
2259	! twice the extinction by the above layers).
2260	!Iabs(iainia) =swdown(iainia) * W_to_mol * RG_to_PAR *&
2261	! Light_Abs_Tot(iainia,jv,ilev) * &
2262	! light_tran_to_level(iainia,jv,ilev) / &
2263	! lai_per_level(iainia,jv,ilev)
2264	Iabs(iainia) =swdown(iainia) * W_to_mol * RG_to_PAR *&
2265	Light_Abs_Tot(iainia,jv,ilev) / &
2266	lai_per_level(iainia,jv,ilev)
2267
2268
2269	! If we have a very dense canopy, the lower levels might
2270	! get almost no light. If the light is close to zero, it
2271	! causes a numerical problem below, as x1 will be zero. Of
2272	! course, if the absorbed light is really small
2273	! photosynthesis will not happen, so we can just skip this
2274	! level. This number is fairly arbitrary, but it should be
2275	! as small as possible.
2276	IF(Iabs(iainia) .LT. 1e-12)CYCLE
2277
2278	! eq. 4 of Yin et al (2009)
2279	! Note that Iabs is very sensitive to the number of layers
2280	! use in the photosynthesis module and to the canopy
2281	! structure. The first sensitivity is related to
2282	! none-linearities to the light conditions and is an
2283	! numerical issue. The sensitivity to the canopy structure
2284	! is a wanted behavior (that's why canopy structure was
2285	! implemented) but makes a substantial difference for the
2286	! total GPP in tropical compared to temperate forests. An
2287	! easy way to prescribe a different canopy structure is by
2288	! changing the range of the prescribed height. The parameter
2289	! ::alpha_LL describes the relationship between Jmax and
2290	! light saturation. It is thought to depend on the amount of
2291	! chlorophyl in the leaves which in turn depends on the
2292	! N-content. For the moment we use a prescribed value. After
2293	! the introduction of the structured canopy we noticed an
2294	! important sensitivity of GPP to this parameter which may
2295	! indicate that a substantial part of our canopy is under
2296	! light saturation.
2297	Jmax(iainia) = vj2(iainia,ilev)
2298	JJ(iainia,ilev) = ( alpha_LL(jv) * Iabs(iainia) + Jmax(iainia) - &
2299	sqrt((alpha_LL(jv) * Iabs(iainia) + Jmax(iainia) )**2. &
2300	- 4 * theta(jv) * Jmax(iainia) * alpha_LL(jv) * &
2301	Iabs(iainia)) ) &
2302	/ ( 2 * theta(jv))
2303
2304	! Store JJ for recalculation in mleb_diffuco_co2. To avoid expensive
2305	! calculations with analytic solutions and loops, JJ will be out.
2306	JJ_out(iainia,jv,ilev) = JJ(iainia,ilev)
2307
2308	IF (printlev_loc>=4 .AND. test_pft == jv) THEN
2309	WRITE(numout,) '************************'
2310	WRITE(numout,*) 'jv,ilev: ',jv,ilev
2311	WRITE(numout,*) 'JJ(iainia,ilev),Jmax(iainia),Iabs(iainia),Rd(iainia,ilev): ',&
2312	JJ(iainia,ilev),Jmax(iainia),Iabs(iainia),Rd(iainia,ilev)
2313	WRITE(numout,*) 'vj2(iainia,ilev),vc2(iainia,ilev), N_profile:',&
2314	vj2(iainia,ilev),vc2(iainia,ilev),N_profile
2315	WRITE(numout,*) 'Light_Abs_Tot(iainia,jv,ilev),Light_Tran_Tot(iainia,jv,ilev):',&
2316	Light_Abs_Tot(iainia,jv,ilev), &
2317	Light_Tran_Tot(iainia,jv,ilev)
2318
2319	WRITE(numout,*) 'water_lim(iainia),T_Rd(iainia):',&
2320	water_lim(iainia),T_Rd(iainia)
2321	WRITE(numout,*) 'vc(iainia),vj(iainia), vcmax(iainia,jv):',&
2322	vc(iainia),vj(iainia),vcmax(iainia,jv)
2323	WRITE(numout,) '************************'
2324	ENDIF
2325
2326	IF ( is_c4(jv) ) THEN
2327	!
2328	! @addtogroup Photosynthesis
2329	! @{
2330	!
2331	! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n
2332	!! \latexonly
2333	!! \input{diffuco_trans_co2_2.4.2.tex}
2334	!! \endlatexonly
2335	! @}
2336	!
2337	! Analytical resolution of the Assimilation based
2338	! Yin et al. (2009)
2339	! Eq. 28 of Yin et al. (2009)
2340	fcyc= 1. - ( 4.(1.-fpsir(jv))(1.+fQ(jv)) + &
2341	3.h_protons(jv)fpseudo(jv) ) / &
2342	( 3.h_protons(jv) - 4.(1.-fpsir(jv)))
2343
2344	! See paragraph after eq. (20b) of Yin et al.
2345	Rm=Rd(iainia,ilev)/2.
2346
2347	! We assume that cs_star equals ci_star
2348	! (see Comment in legend of Fig. 6 of Yin et al. (2009)
2349	! Equation 26 of Yin et al. (2009)
2350	Cs_star = (gbs(jv) * low_gamma_star(iainia) * &
2351	Oi - ( 1. + low_gamma_star(iainia) * alpha(jv) / &
2352	0.047) * Rd(iainia,ilev) + Rm ) / ( gbs(jv) + kp(jv) )
2353
2354	! eq. 11 of Yin et al (2009)
2355	J2 = JJ(iainia,ilev) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) )
2356
2357	! Equation right after eq. (20d) of Yin et al. (2009)
2358	z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc ))
2359
2360	VpJ2 = fpsir(jv) * J2 * z / 2.
2361
2362	A_3=9999.
2363
2364	! See eq. right after eq. 18 of Yin et al. (2009)
2365	! The loop over limit_photo calculates the assimultion
2366	! rate of the Rubisco-limited CO2 assimilation (Ac) and
2367	! the rate of electron transport-limited CO2 assimultion
2368	! (Aj). The overall assimilation will be the slower
2369	! of these two processes.
2370	DO limit_photo=1,2
2371	! Is Vc limiting the Assimilation
2372	IF ( limit_photo .EQ. 1 ) THEN
2373	a = 1. + kp(jv) / gbs(jv)
2374	b = 0.
2375	x1 = vc2(iainia,ilev)
2376	x2 = KmC(iainia)/KmO(iainia)
2377	x3 = KmC(iainia)
2378	! Is J limiting the Assimilation
2379	ELSE
2380	a = 1.
2381	b = VpJ2
2382	x1 = (1.- fpsir(jv)) * J2 * z / 3.
2383	x2 = 7. * low_gamma_star(iainia) / 3.
2384	x3 = 0.
2385	ENDIF
2386
2387	m=fvpd(iainia)-g0var(iainia)/gb_co2(iainia)
2388	d=g0var(iainia)(Ca(iainia)-Cs_star) + fvpd(iainia)Rd(iainia,ilev)
2389	f=(b-Rm-low_gamma_star(iainia)Oigbs(jv))x1d + agbs(jv)x1Ca(iainia)d
2390	j=(b-Rm+gbs(jv)x3 + x2gbs(jv)Oi)m + (alpha(jv)x2/0.047-1.)d &
2391	+ agbs(jv)(Ca(iainia)*m - d/gb_co2(iainia) - (Ca(iainia) - Cs_star ))
2392
2393	g=(b-Rm-low_gamma_star(iainia)Oigbs(jv))x1m - (alpha(jv)low_gamma_star(iainia)/0.047+1.)x1*d &
2394	+ agbs(jv)x1(Ca(iainia)m - d/gb_co2(iainia) - (Ca(iainia)-Cs_star ))
2395
2396	h=-((alpha(jv)low_gamma_star(iainia)/0.047+1.)x1m + (agbs(jv)x1(m-1.))/gb_co2(iainia) )
2397	i= ( b-Rm + gbs(jv)x3 + x2gbs(jv)Oi )d + agbs(jv)Ca(iainia)*d
2398	l= ( alpha(jv)x2/0.047 - 1.)m - (agbs(jv)(m-1.))/gb_co2(iainia)
2399
2400	p = (j-(h-l*Rd(iainia,ilev))) / l
2401	q = (i+j*Rd(iainia,ilev)-g) / l
2402	r = -(f-i*Rd(iainia,ilev)) / l
2403
2404	! See Yin et al. (2009) and Baldocchi (1994)
2405	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
2406	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
2407
2408
2409	IF (QQ .GE. 0._r_std) THEN
2410	IF (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) THEN
2411	PSI = ACOS(UU/(QQ**1.5_r_std))
2412	A_3_tmp = -2._r_std * SQRT(QQ) * COS(( PSI + 4._r_std * PI)/3._r_std ) - p / 3._r_std
2413	IF (( A_3_tmp .LT. A_3 )) THEN
2414	A_3 = A_3_tmp
2415	info_limitphoto(iainia,jv,ilev)=2.
2416	ELSE
2417	! In case, J is not limiting the assimilation
2418	! we have to re-initialise a, b, x1, x2 and x3 values
2419	! in agreement with a Vc-limited assimilation
2420	a = 1. + kp(jv) / gbs(jv)
2421	b = 0.
2422	x1 = vc2(iainia,ilev)
2423	x2 = KmC(iainia)/KmO(iainia)
2424	x3 = KmC(iainia)
2425	info_limitphoto(iainia,jv,ilev)=1.
2426	ENDIF
2427	ENDIF
2428	ELSE
2429	WRITE(numout,*) 'WARNING: unexpected value for QQ (diffuco.f90)'
2430	ENDIF
2431
2432	IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(iainia,ilev)) ) ) THEN
2433	IF ( printlev>=4 ) THEN
2434	WRITE(numout,*) 'We have a problem in diffuco_trans_co2'
2435	WRITE(numout,*) 'no real positive solution found for pft:',jv
2436	WRITE(numout,*) 'temp_air:',temp_air(iainia)
2437	WRITE(numout,*) 'vpd:',vpd(iainia)
2438	WRITE(numout,*) 'Rd:',Rd(iainia,ilev)
2439	END IF
2440	A_3 = -Rd(iainia,ilev)
2441	ENDIF
2442
2443	assimi(iainia,jv,ilev) = A_3
2444
2445	IF(printlev_loc>=4 .AND. jv == test_pft .AND. iainia == test_pft)THEN
2446	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2447	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2448	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2449	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2450	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2451	ENDIF
2452
2453	!++++ TEMP +++++
2454	IF(( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) )) == zero)THEN
2455	! This causes a problem, since we really should not
2456	! be dividing by zero. However, it seems to happen
2457	! fairly rarely, so all we'll do is cycle to the
2458	! next loop and keep track of this warning so we
2459	! can look at it afterwards.
2460	warnings(iainia,jv,iwphoto)=warnings(iainia,jv,iwphoto)+un
2461
2462	IF(err_act.GT.1)THEN
2463	WRITE(numout,*) 'Diffuco, Getting ready to divide by zero! C4'
2464	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2465	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2466	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2467	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2468	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2469	ENDIF
2470
2471	CYCLE
2472
2473	ENDIF
2474
2475	IF ( ABS( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) .LT. min_sechiba ) THEN
2476	gs(iainia,ilev) = g0var(iainia)
2477	!leaf_ci keeps its initial value (Ca)
2478	ELSE
2479	! Eq. 24 of Yin et al. (2009)
2480	Obs = ( alpha(jv) * assimi(iainia,jv,ilev) ) / ( 0.047 * gbs(jv) ) + Oi
2481	! Eq. 23 of Yin et al. (2009)
2482	Cc_loc(iainia,ilev) = ( ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) * &
2483	( x2 * Obs + x3 ) + low_gamma_star(iainia) * Obs * x1 ) &
2484	/ MAX(min_sechiba, x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2485	! Eq. 22 of Yin et al. (2009)
2486	leaf_ci(iainia,jv,ilev) = ( Cc_loc(iainia,ilev) - ( b - assimi(iainia,jv,ilev) - Rm ) / gbs(jv) ) / a
2487
2488	!This is eq. 25 of Yin et al. (2009), but we want to
2489	!use Ca instead of Cs. Therefore we need to take into
2490	!account the CO2 transfer along the path from Ca to Cs,
2491	!thus accounting for the boundary layer conductance (gs).
2492	gs(iainia,ilev) = g0var(iainia) +(assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) / &
2493	((Ca(iainia)-assimi(iainia,jv,ilev))/gb_co2(iainia)-Cs_star ) * fvpd(iainia)
2494	ENDIF
2495	!
2496	ENDDO !! ENDDO about limit photo
2497	!
2498	ELSE
2499	!
2500	! @addtogroup Photosynthesis
2501	! @{
2502	!
2503	! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
2504	!! \latexonly
2505	!! \input{diffuco_trans_co2_2.4.3.tex}
2506	!! \endlatexonly
2507	! @}
2508	!
2509	!
2510	!
2511	A_1=9999.
2512
2513	! See eq. right after eq. 18 of Yin et al. (2009)
2514	DO limit_photo=1,2
2515	! Is Vc limiting the Assimilation
2516	IF ( limit_photo .EQ. 1 ) THEN
2517
2518	IF (printlev_loc>=4) THEN
2519	WRITE(numout,*) 'IF - A_1, Do in limit_photo limit_photo .EG. 1 '
2520	ENDIF
2521	x1 = vc2(iainia,ilev)
2522	! It should be O not Oi (comment from Vuichard)
2523	x2 = KmC(iainia) * ( 1. + 2gamma_star(iainia)Sco(iainia) / KmO(iainia) )
2524	! Is J limiting the Assimilation
2525	ELSE
2526	x1 = JJ(iainia,ilev)/4.
2527	x2 = 2. * gamma_star(iainia)
2528	ENDIF
2529
2530
2531	! See Appendix B of Yin et al. (2009)
2532	a = g0var(iainia) * ( x2 + gamma_star(iainia) ) + &
2533	( g0var(iainia) / gm(iainia) + fvpd(iainia) ) * ( x1 - Rd(iainia,ilev) )
2534	b = Ca(iainia) * ( x1 - Rd(iainia,ilev) ) - gamma_star(iainia) * x1 - Rd(iainia,ilev) * x2
2535	c = Ca(iainia) + x2 + ( 1./gm(iainia) + 1./gb_co2(iainia) ) * ( x1 - Rd(iainia,ilev) )
2536	d = x2 + gamma_star(iainia) + ( x1 - Rd(iainia,ilev) ) / gm(iainia)
2537	m = 1./gm(iainia) + ( g0var(iainia)/gm(iainia) + fvpd(iainia) ) * ( 1./gm(iainia) + 1./gb_co2(iainia) )
2538
2539	p = - ( d + (x1 - Rd(iainia,ilev) ) / gm(iainia) + a * (1./gm(iainia) + 1./gb_co2(iainia) ) + &
2540	( g0var(iainia)/gm(iainia) + fvpd(iainia)) * c) / m
2541
2542	q = ( d * ( x1 - Rd(iainia,ilev) ) + ac + ( g0var(iainia)/gm(iainia) + fvpd(iainia) ) b ) / m
2543	r = - a * b / m
2544
2545	! See Yin et al. (2009)
2546	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
2547	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
2548
2549	IF (QQ .GE. 0._r_std) THEN
2550
2551	IF (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) THEN
2552	PSI = ACOS(UU/(QQ**1.5_r_std))
2553	A_1_tmp = -2._r_std * SQRT(QQ) * COS( PSI / 3._r_std ) - p / 3._r_std
2554
2555	IF (( A_1_tmp .LT. A_1 )) THEN
2556	A_1 = A_1_tmp
2557	info_limitphoto(iainia,jv,ilev)=2.
2558	ELSE
2559	! In case, J is not limiting the assimilation
2560	! we have to re-initialise x1 and x2 values
2561	! in agreement with a Vc-limited assimilation
2562	x1 = vc2(iainia,ilev)
2563	! It should be O not Oi (comment from Vuichard)
2564	x2 = KmC(iainia) * ( 1. + 2gamma_star(iainia)Sco(iainia) / KmO(iainia) )
2565	info_limitphoto(iainia,jv,ilev)=1.
2566	ENDIF
2567	ENDIF
2568	ENDIF
2569	ENDDO
2570
2571	IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(iainia,ilev)) ) ) THEN
2572	IF ( printlev>=4 ) THEN
2573	WRITE(numout,*) 'We have a problem in diffuco_trans_co2'
2574	WRITE(numout,*) 'no real positive solution found for pft:',jv
2575	WRITE(numout,*) 'temp_air:',temp_air(iainia)
2576	WRITE(numout,*) 'vpd:',vpd(iainia)
2577	WRITE(numout,*) 'Setting the solution to -Rd(iainia,ilev):',-Rd(iainia,ilev)
2578	END IF
2579	A_1 = -Rd(iainia,ilev)
2580	ENDIF
2581	assimi(iainia,jv,ilev) = A_1
2582
2583	! Eq. 18 of Yin et al. (2009)
2584
2585	IF(printlev_loc>=4 .AND. jv == test_pft .AND. iainia == test_pft)THEN
2586	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2587	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2588	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2589	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2590	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2591	ENDIF
2592
2593	!++++ TEMP +++++
2594	IF(( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) )) == zero)THEN
2595
2596	! This causes a problem, since we really should not
2597	! be dividing by zero. However, it seems to happen
2598	! fairly rarely, so all we'll do is cycle to the
2599	! next loop and keep track of this warning so we
2600	! can look at it afterwards.
2601	warnings(iainia,jv,iwphoto)=warnings(iainia,jv,iwphoto)+un
2602
2603	IF(err_act.GT.1)THEN
2604	WRITE(numout,*) 'Diffuco, Getting ready to divide by zero! C3'
2605	WRITE(numout,*) '( x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) ', &
2606	( x1 - (assimi(iainia,jv,ilev) + Rd(iainia,ilev) ))
2607	WRITE(numout,*) 'x1,assimi(iainia,jv,ilev),Rd(iainia,ilev): ', &
2608	x1,assimi(iainia,jv,ilev),Rd(iainia,ilev)
2609	WRITE(numout,*) 'jv,ilev,iainia: ', jv,ilev,iainia
2610	ENDIF
2611
2612	CYCLE
2613
2614	ENDIF
2615	!++++++++++++++++++++++++++++++
2616
2617	IF ( ABS( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) .LT. min_sechiba ) THEN
2618	gs(iainia,ilev) = g0var(iainia)
2619	ELSE
2620	! Eq. 18 of Yin et al. (2009)
2621	Cc_loc(iainia,ilev) = ( gamma_star(iainia) * x1 + ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) * x2 ) &
2622	/ MAX( min_sechiba, x1 - ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) )
2623	! Eq. 17 of Yin et al. (2009)
2624	leaf_ci(iainia,jv,ilev) = Cc_loc(iainia,ilev) + assimi(iainia,jv,ilev) / gm(iainia)
2625	! See eq. right after eq. 15 of Yin et al. (2009)
2626	ci_star = gamma_star(iainia) - Rd(iainia,ilev) / gm(iainia)
2627	!
2628	! Eq. 15 of Yin et al. (2009)
2629	gs(iainia,ilev) = g0var(iainia) + ( assimi(iainia,jv,ilev) + Rd(iainia,ilev) ) / ( leaf_ci(iainia,jv,ilev) &
2630	- ci_star ) * fvpd(iainia)
2631	ENDIF
2632
2633
2634	ENDIF !C3 vs C4
2635
2636	!
2637	!
2638	! this is only for output to netcdf
2639	gs_diffuco_output(iainia,jv,ilev) = gs(iainia,ilev)
2640	!
2641	!
2642	! @addtogroup Photosynthesis
2643	! @{
2644	!
2645	!! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n
2646	!! \latexonly
2647	!! \input{diffuco_trans_co2_2.4.4.tex}
2648	!! \endlatexonly
2649	! @}
2650	!
2651	!
2652	! keep stomatal conductance of the levels which we decide
2653	! are the top levels. Notice that we multiply by the LAI
2654	! in the level here, so this is no longer a per leaf quantity,
2655	! but the whole quantity.
2656	IF ( top_level(iainia,ilev) ) THEN
2657	gs_top(iainia) = gs_top(iainia)+gs(iainia,ilev)*lai_per_level(iainia,jv,ilev)
2658	!
2659	ENDIF
2660	!
2661	! @addtogroup Photosynthesis
2662	! @{
2663	!
2664	!! 2.4.5 Integration at the canopy level\n
2665	!! \latexonly
2666	!! \input{diffuco_trans_co2_2.4.5.tex}
2667	!! \endlatexonly
2668	! @}
2669	! total assimilation and conductance
2670	!
2671	!++++++++++++++++++++
2672	! The equations give the amount of carbon produced per m**2 if
2673	! leaf surface. Therefore, since we want GPP in units of m**2 of land surface,
2674	! we have to multiply by the LAI.
2675
2676	assimtot(iainia,jv) = assimtot(iainia,jv) + &
2677	assimi(iainia,jv,ilev) * lai_per_level(iainia,jv,ilev)
2678
2679	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
2680	profile_assimtot(iainia,ilev) = profile_assimtot(iainia,ilev) + &
2681	& assimi(iainia,jv,ilev) * lai_per_level(iainia,jv,ilev)
2682	profile_Rdtot(iainia,ilev) = profile_Rdtot(iainia,ilev) + &
2683	& Rd(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2684	profile_gstot(iainia,ilev) = profile_gstot(iainia,ilev) + &
2685	& gs(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2686	END IF ! (ok_impose_can_structure ) THEN
2687
2688	! assimi(iainia)
2689
2690	Rdtot(iainia) = Rdtot(iainia) + &
2691	Rd(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2692	gstot(iainia) = gstot(iainia) + &
2693	gs(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2694	assimi_lev(iainia,jv,ilev) = assimi(iainia,jv,ilev)
2695
2696	IF(test_pft == jv .AND. printlev_loc>=4)THEN
2697	WRITE(numout,'(A,2I5,10ES20.8)') 'assimi middle',ilev,&
2698	iainia,assimi(iainia,jv,ilev),assimtot(iainia,jv),gs(iainia,ilev),&
2699	gstot(iainia),lai_per_level(iainia,jv,ilev)
2700	ENDIF
2701
2702	!! (JR 100714) we need to retain the gs (stomatal conductance) per level and per PFT for use
2703	!! outside of this routine, in the case that the supply term constraint is applied and the
2704	!! stomatal conductance profile needs to be re-calculated.
2705	gs_store(iainia, jv, ilev) = gs(iainia,ilev)
2706	!
2707
2708	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
2709	profile_gstot(iainia, ilev) = profile_gstot(iainia, ilev) + &
2710	& gs(iainia,ilev) * lai_per_level(iainia,jv,ilev)
2711	END IF ! (ok_impose_can_structure ) THEN
2712
2713	!! (JR 100714) we also retain the gstot (stomatal conductance over all canopy levels, constrained
2714	!! by LAI, again for the recalculation of the stomatal conductance profile.
2715	gstot_store(iainia, jv) = gstot_store(iainia, jv) + (gs(iainia,ilev) * lai_per_level(iainia,jv,ilev))
2716	!
2717	ENDIF ! if there is absorbed light
2718
2719	ENDDO ! loop over nlevels_tot
2720
2721	!! Calculated intercellular CO2 over nlai needed for the chemistry module
2722	cim(iainia,jv) = zero
2723	laisum(iainia) = zero
2724	DO ilev = 1,nlevels_tot
2725	IF(lai_per_level(iainia,jv,ilev).GT.zero .AND. leaf_ci(iainia,jv,ilev).GT.zero)THEN
2726	cim(iainia,jv)= cim(iainia,jv)+leaf_ci(iainia,jv,ilev)*lai_per_level(iainia,jv,ilev)
2727	laisum(iainia)=laisum(iainia) + lai_per_level(iainia,jv,ilev)
2728	ENDIF
2729	ENDDO
2730	IF (laisum(iainia).GT.zero) THEN
2731	cim(iainia,jv)= cim(iainia,jv)/laisum(iainia)
2732	ELSE
2733	cim(iainia,jv)=zero
2734	ENDIF
2735
2736	ENDDO assim_loop ! loop over assimulation points
2737	IF(jv==testpft) THEN
2738	! CALL histwrite_p(hist_id, 'Cc_loc', kjit, Cc_loc, kjpindex*(nlevels_tot+1), indexlai)
2739	! CALL histwrite_p(hist_id, 'Vc', kjit, Vc2, kjpindex*(nlevels_tot+1), indexlai)
2740	! CALL histwrite_p(hist_id, 'Vj', kjit, JJ, kjpindex*(nlevels_tot+1), indexlai)
2741	! CALL histwrite_p(hist_id, 'limitphoto', kjit, info_limitphoto, kjpindexnvm(nlevels_tot+1), indexlai)
2742	! CALL histwrite_p(hist_id, 'gammastar', kjit, gamma_star, kjpindex, index)
2743	! CALL histwrite_p(hist_id, 'Kmo', kjit, Kmo, kjpindex, index)
2744	! CALL histwrite_p(hist_id, 'Kmc', kjit, Kmc, kjpindex, index)
2745	! CALL histwrite_p(hist_id, 'gm', kjit, gm, kjpindex, index)
2746	! CALL histwrite_p(hist_id, 'gs', kjit, gs, kjpindex*(nlevels_tot+1), indexlai)
2747	! CALL histwrite_p(hist_id, 'assimi', kjit, assimi, kjpindexnvm(nlevels_tot+1), indexlai)
2748	! CALL histwrite_p(hist_id, 'Rd', kjit, Rd, kjpindex*(nlevels_tot+1), indexlai)
2749	ENDIF
2750
2751	!
2752	!! 2.5 Calculate resistances
2753
2754	DO inia=1,nia
2755
2756	!
2757	iainia=index_assi(inia)
2758	!
2759
2760	! We hit a fringe case above where we could not
2761	! calculate the GPP. Skip this loop so we don't
2762	! divide by zero.
2763	IF(lskip(iainia)) CYCLE
2764
2765
2766	!! Mean stomatal conductance for CO2 (mol m-2 s-1)
2767	gsmean(iainia,jv) = gstot(iainia)
2768	!
2769	! cimean is the "mean ci" calculated in such a way that assimilation
2770	! calculated in enerbil is equivalent to assimtot
2771	!
2772
2773	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
2774	DO ilev = 1, nlevels_tot
2775	profile_gsmean(iainia, jv, ilev) = profile_gstot(iainia, ilev)
2776	IF(printlev_loc >= 4) WRITE(numout,*) '110215 profile_gsmean by level: ', &
2777	ilev, ' is: ', profile_gsmean(iainia, jv, ilev)
2778	ENDDO ! ilev = 1, nlevels_tot
2779	END IF ! (ok_impose_can_structure ) THEN
2780
2781
2782	! cimean is completely diagnostic and not
2783	! used for anything. If the analytical photosynthesis is unable to find
2784	! a solution for the A_1 coefficients in any of the levels, assimtot will
2785	! be equal to -Rd, which means we will be dividing by zero. This loop
2786	! is to prevent that, and since cimean is diagnostic is doesn't matter
2787	! what the value is.
2788	IF ( ABS(gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) .GT. min_sechiba) THEN
2789	cimean(iainia,jv) = (fvpd(iainia)*(assimtot(iainia,jv)+Rdtot(iainia))) /&
2790	(gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) + gamma_star(iainia)
2791	ELSE
2792	cimean(iainia,jv) = gamma_star(iainia)
2793	ENDIF
2794
2795
2796	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
2797	DO ilev = 1, nlevels_tot
2798	IF((profile_assimtot(iainia,ilev)+profile_Rdtot(iainia,ilev)) .GT. min_sechiba)THEN
2799	profile_cimean(iainia,jv,ilev) = (profile_gs_store(iainia,jv,ilev)-g0(jv)) &
2800	/ (fvpd(iainia)*(profile_assimtot(iainia,ilev)+profile_Rdtot(iainia,ilev))) &
2801	+ gamma_star(iainia)
2802	IF(printlev_loc >= 4) THEN
2803	WRITE(numout,*) '210515 fvpd:', fvpd(iainia)
2804	WRITE(numout,*) 'profile_assimtot:', profile_assimtot(iainia,ilev)
2805	WRITE(numout,*) 'profile Rdtot:', profile_Rdtot(iainia,ilev)
2806	ENDIF ! IF(printlev >= 4)
2807	ELSE
2808	! According to Nicolas, this variable is completely diagnostic and not
2809	! used for anything. If the analytical photosynthesis is unable to find
2810	! a solution for the A_1 coefficients in any of the levels, assimtot will
2811	! be equal to -Rd, which means we will be dividing by zero. This loop
2812	! is to prevent that, and since cimean is diagnostic is doesn't matter
2813	! what the value is.
2814	profile_cimean(iainia,jv,ilev) = val_exp
2815	ENDIF
2816	! WRITE(numout, *) '100215b profile_cimean by level: ', ilev, ' is: ', profile_cimean(iainia, jv, ilev)
2817	ENDDO ! ilev = 1, nlevels_tot
2818	END IF ! (ok_impose_can_structure ) THEN
2819
2820
2821	! conversion from umol m-2 (PFT) s-1 to gC m-2 (mesh area) tstep-1
2822	gpp(iainia,jv) = assimtot(iainia,jv)12e-6veget_max(iainia,jv)*dt_sechiba
2823	istep=istep+1
2824
2825	IF (printlev_loc>=4) THEN
2826	WRITE(numout,*) 'GPP DIFFUCO: ',istep,gpp(iainia,jv),iainia,jv,assimtot(iainia,jv),veget_max(iainia,jv),dt_sechiba
2827	ENDIF
2828
2829	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
2830	DO ilev = 1, nlevels_tot
2831	profile_gpp(iainia,jv,ilev) = profile_assimtot(iainia, ilev)12e-6veget_max(iainia,jv)*dt_sechiba
2832	IF(printlev_loc >= 4) WRITE(numout, *) '100215b profile_gpp by level: ', &
2833	ilev, ' is: ', profile_gpp(iainia, jv, ilev)
2834	END DO ! ilev = 1, nlevels_tot
2835	END IF ! (ok_impose_can_structure ) THEN
2836
2837	! As in Pearcy, Schulze and Zimmermann
2838	! Measurement of transpiration and leaf conductance
2839	! Chapter 8 of Plant Physiological Ecology
2840	! Field methods and instrumentation, 1991
2841	! Editors:
2842	!
2843	! Robert W. Pearcy,
2844	! James R. Ehleringer,
2845	! Harold A. Mooney,
2846	! Philip W. Rundel
2847	!
2848	! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online)
2849	!
2850	!
2851	!
2852	!++++++++ CHECK +++++++++
2853	! Do we need to normalise gs_top here by the LAI in the level?
2854	! In the trunk, it seems like they use an LAI of 0.1, so renormalise
2855	! to that. The default value of lai_top is 0.1.
2856	!
2857	!! (JR 100714)
2858	!! gs_top is normalised in the trunk over the entire canopy LAI, as it
2859	!! meant to provide backwards compatability with the 'old' method for calculating
2860	!! the stomatal conductance (in diffuco_trans) that calculates gs according to
2861	!! incident SW radiation (rather than CO2 concentration, as here in diffuco_trans_co2)
2862	!!
2863	!! rstruct in diffuco_trans is an extra resistance term added to reflect the fact
2864	!! that less light reaches the lower canopy, and the total resistance used to calculate
2865	!! vbeta3, and hence the transpiration, is (rstruct + rveget)
2866	!!
2867	!! The rstruct term makes the code more difficult to follow, and is no longer has
2868	!! a physical requirement. We are getting small values of gstop as calculated below,
2869	!! which was resulting in big resistances when rstruct = min_sechiba and
2870	!! rveget = un/gstop and total resistance = (rstruct + rveget) - see commented out
2871	!! code below.
2872	!!
2873	!! The most transparent structure is to set total resistance = rveget = (1 / gstot)
2874	!! This is analogous to the way in which assimtot, and hence canopy GPP, is calculated
2875	!! above
2876	gs_top(iainia)=gs_top(iainia)/totlai_top*lai_top(jv)
2877	!+++++++++++++++++++++++++
2878
2879	gstot(iainia) = mol_to_m_1 (temp_air(iainia)/tp_00)&
2880	(pb_std/pb(iainia))gstot(iainia)ratio_H2O_to_CO2
2881	gstop(iainia) = mol_to_m_1 * (temp_air(iainia)/tp_00)*&
2882	(pb_std/pb(iainia))gs_top(iainia)ratio_H2O_to_CO2
2883
2884
2885	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
2886	DO ilev = 1, nlevels_tot
2887
2888	profile_gsmean(iainia,jv,ilev) = profile_gsmean(iainia,jv,ilev) * 1e-6
2889
2890	profile_gstot(iainia,ilev) = mol_to_m_1 (temp_air(iainia)/tp_00)&
2891	(pb_std/pb(iainia))profile_gstot(iainia,ilev)ratio_H2O_to_CO2
2892
2893	profile_gstop(iainia,ilev) = mol_to_m_1 * (temp_air(iainia)/tp_00)*&
2894	(pb_std/pb(iainia))gs_top(iainia)ratio_H2O_to_CO2*&
2895	lai_per_level(iainia,jv,ilev)
2896
2897	END DO ! ilev = 1, nlevels_tot
2898	END IF ! (ok_impose_can_structure ) THEN
2899
2900
2901	!IF (jv == test_pft .AND. test_grid == 1) THEN
2902	! WRITE(numout, '(A20)' ) 'Without water stress:'
2903	! WRITE(numout, '(A20, I5)' ) 'pft type is: ', jc
2904	! WRITE(numout, '(A20, ES15.5)') 'Old_Gs_total:', gstot(iainia)
2905	!ENDIF
2906
2907	! THIS PART OF CORRECTION OF ::gstot SHOULD BE CONSISTENCY WITH
2908	! "SECHIBA_HYDROL_ARCH" where the ::gstot is recalculated, again.
2909
2910	!++++++ CHECK ++++++
2911	! We introduce ::coupled_lai to calculate the :: vbeta3
2912	! In an closed canopy not all the leaves fully interact with the
2913	! atmosphere because leaves can shelter each other. In a more open
2914	! canopy most leaves can interact with the atmosphere. We use the
2915	! ratio of lai over ::Pgap as a proxy for the amount of canopy that
2916	! can interact with the atmosphere. Clearly that amount of lai that
2917	! interacts should never exceed the:: lai.
2918
2919	! Use the lai also as coupled lai
2920	coupled_lai(iainia, jv) = lai(iainia, jv)
2921
2922	!---TEMP---
2923	!!$ IF (printlev_loc>=4) THEN
2924	!!$ WRITE(numout,*) 'coupled LAI, ', jv, coupled_lai(iainia, jv)
2925	!!$ WRITE(numout,*) 'total LAI, ', jv, lai(iainia, jv)
2926	!!$ WRITE(numout,*) 'veget, ', jv, veget(iainia, jv)
2927	!!$ WRITE(numout,*) 'biomass(ileaf), ',biomass(iainia, jv, ileaf, icarbon)
2928	!!$ ENDIF
2929	!----------
2930
2931	!+++++CHECK++++
2932	!Here, we try to use coupled_lai to reset/refine the total stomata conductance
2933	!which interacted with the atmosphere for the transpiration for each pft.
2934	IF (lai(iainia,jv) > min_sechiba) THEN
2935	gstot(iainia) = gstot(iainia)* (coupled_lai(iainia,jv)/lai(iainia,jv))
2936	ELSE
2937	gstot(iainia) = 0
2938	END IF
2939	!++++++++++++++
2940
2941	!+++++++++CHECK++++++++
2942	! To get the resonable transpiration we had to set :: lai_top to
2943	! ::lai
2944	! so we were using the entire canopy. Therefore we no longer use
2945	! ::lai_top and simply make use of the entire canopy to calculate the
2946	! vegetation resistance from inverse of ::gstot so, we calculate the
2947	! :: rveget using :: 1/gstot. This case where gstot(iana)=zero hints
2948	! at another problem: we have vegetation but no photosynthesis.
2949	IF (gstot(iainia).GT.zero)THEN
2950	rveget(iainia, jv) = un / gstot(iainia)
2951	ELSE
2952	rveget(iainia, jv) = un
2953	CALL ipslerr(2,'diffuco_trans_co2','gstot is zero! This hints at another problem',&
2954	'do we have canopy layers with LAI but without absorbed light?','')
2955	ENDIF
2956	!+++++++++++++++++++++++++++++++++++++++
2957
2958	!rveget(iainia,jv) = un/gstop(iainia)
2959	!!$ rveget_min(iainia,jv) = (defc_plus / kzero(jv)) * &
2960	!!$ (un / SUM(lai_per_level(iainia,jv,:)))
2961	!
2962	! rstruct is the difference between rtot (=1./gstot) and rveget
2963	!
2964	! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
2965	!rstruct(iainia,jv) = un/gstot(iainia) - &
2966	! rveget(iainia,jv)
2967
2968	!++++++++ CHECK +++++++
2969	!It seems that the calcualtion of rstruct is duplicated with the calculation of rvegt
2970	!So, we set the ::rstruct as zero
2971	!rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
2972	! rveget(iainia,jv), min_sechiba)
2973	!
2974	!
2975	!! wind is a global variable of the diffuco module.
2976	speed = MAX(min_wind, wind(iainia))
2977	!
2978	! beta for transpiration
2979	!
2980	! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
2981	!! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
2982	!! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
2983	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2984	! (un - zqsvegrap(iainia)) * &
2985	! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
2986	! rstruct(iainia,jv))))
2987	!! Global resistance of the canopy to evaporation
2988	!cresist=(un / (un + speed * q_cdrag(iainia) * &
2989	! veget(iainia,jv)/veget_max(iainia,jv) * &
2990	! (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
2991
2992	! (JR 100714)
2993	! previously the resistance term in this expression was (rveget + rstruct), but we now
2994	! use rveget alone
2995
2996	cresist(iainia,jv)=(un / (un + speed * q_cdrag(iainia) * &
2997	veget(iainia,jv)/veget_max(iainia,jv) * &
2998	(rveg_pft(jv)*(rveget(iainia,jv)))))
2999
3000
3001	! cresist=un / (un + speed * q_cdrag(iainia) * &
3002	! veget(iainia,jv)/veget_max(iainia,jv) * &
3003	! (rveg_pft(jv)*rveget(iainia,jv)))
3004
3005	!+++CHECK+++
3006	!For the interception we assume -as earlier- that the surface area
3007	!that intercepts is well estimated by ::veget. This is wrong because
3008	!every leav can contribute to the interception. Depending on weather
3009	!the parameters for interception are bulk values or single leave. Our
3010	!approach is correct or needd to be improved.
3011	!
3012	IF ( humrel(iainia,jv) >= min_sechiba ) THEN
3013	vbeta3(iainia,jv) = veget(iainia,jv) * &
3014	(un - zqsvegrap(iainia)) * cresist(iainia,jv) + &
3015	MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
3016	zqsvegrap(iainia) * cresist(iainia,jv) )
3017	ELSE
3018	! Because of a minimum conductance g0, vbeta3 cannot be zero even if humrel=0
3019	! in the above equation.
3020	! Here, we force transpiration to be zero when the soil cannot deliver it
3021	vbeta3(iainia,jv) = zero
3022	END IF
3023
3024	IF (ok_mleb .AND. .NOT. ok_impose_can_structure ) THEN
3025	DO ilev = 1, nlevels_tot
3026
3027	! Fix the gstot equal zero condition to aviod the issue of divided by zero
3028	IF ( profile_gstot(iainia,ilev) .LE. min_sechiba ) THEN
3029	profile_gstot(iainia,ilev) = 1.0d-6
3030	ENDIF
3031
3032	profile_rveget(iainia,jv,ilev) = un/profile_gstot(iainia,ilev)
3033
3034	! need to resolve the rstruct issue
3035	profile_rstruct(iainia,jv,ilev) = MAX (un/profile_gstot(iainia,ilev) - &
3036	& profile_rveget(iainia,jv,ilev), min_sechiba)
3037
3038	profile_speed(ilev) = MAX(min_wind, u_speed(ilev)) ! n.b. need to replace this with another wind profile parameterisation
3039
3040	profile_cresist(ilev) = (un / (un + u_speed(ilev) * q_cdrag(iainia) * &
3041	& (rveg_pft(jv) * (profile_rveget(iainia,jv,ilev) + profile_rstruct(iainia,jv,ilev)) ) ) )
3042
3043	profile_vbeta3(iainia,jv,ilev) = veget_max(iainia,jv) * (un - zqsvegrap(iainia)) * profile_cresist(ilev) + &
3044	& MIN( vbeta23(iainia,jv), &
3045	veget_max(iainia,jv) * zqsvegrap(iainia) * profile_cresist(ilev) )
3046
3047	END DO ! ilev = 1, nlevels_tot
3048	END IF ! (ok_impose_can_structure ) THEN
3049
3050	IF (printlev_loc>=4) THEN
3051	DO ilev = 1, nlevels_tot
3052	WRITE(numout, *) 'profile_rveget(iainia, jv, ilev) is: ', profile_rveget(iainia, jv, ilev)
3053	END DO
3054
3055	DO ilev = 1, nlevels_tot
3056	WRITE(numout, *) 'profile_vbeta3(iainia, jv, ilev) is: ', profile_vbeta3(iainia, jv, ilev)
3057	END DO
3058
3059	DO ilev = 1, nlevels_tot
3060	WRITE(numout, *) 'profile_cresist(ilev) is: ', (profile_cresist(ilev))
3061	END DO
3062	END IF !(printlev_loc)
3063
3064	! vbeta3pot for computation of potential transpiration (needed for irrigation)
3065	vbeta3pot(iainia,jv) = MAX(zero, veget(iainia,jv) * cresist(iainia,jv))
3066	!
3067	!
3068
3069	ENDDO !loop over grid points
3070
3071	!
3072	ENDIF ! if nia LT zero
3073
3074	ENDDO ! loop over vegetation types
3075
3076
3077	DO ic = 1, kjpindex
3078	DO jc = 1, nvm
3079	!++++++TEMP++++++++++++++++++++++++++++++++++++
3080	! IF (printlev_loc>=4 .AND. jc == test_pft) THEN
3081	! WRITE(714, '(A20, ES20.4)') 'Rcanopy,', (un / (un + speed q_cdrag(ic) &
3082	! (rveg_pft(jc)*(rveget(ic,jc) + rstruct(ic,jc)))))
3083
3084	! WRITE(numout,'(A20, ES20.4)') 'Rcanopy',(un / (un + speed q_cdrag(ic) &
3085	! (rveg_pft(jc)*(rveget(ic,jc) + rstruct(ic,jc)))))
3086	! ENDIF
3087	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3088
3089
3090
3091	IF(veget_max(ic,jc) == zero)THEN
3092	! this vegetation type is not present, so no reason to do the
3093	! calculation
3094	CYCLE
3095	ENDIF
3096
3097	DO kc = 1, nlevels_tot
3098	gs_column_total(ic,jc) = gs_column_total(ic,jc) + gs_store(ic,jc,kc) * lai_per_level(ic,jc,kc)
3099
3100	gstot_component(ic,jc,kc) = gs_store(ic,jc,kc) * lai_per_level(ic,jc,kc)
3101
3102	IF (gstot_store(ic,jc) .eq. 0.0d0) THEN
3103	gstot_frac(ic,jc,kc) = 0.0d0
3104	ELSE
3105	gstot_frac(ic,jc,kc) = gstot_component(ic,jc,kc) / gstot_store(ic,jc)
3106	END IF
3107	END DO ! kc = 1, nlevels_tot
3108
3109	END DO ! jc = 1, nvm
3110	END DO ! ic = 1, kjpindex
3111
3112
3113
3114	Abs_light(:,:,:) = zero
3115	Abs_light_can(:,:) = zero
3116
3117	DO ic = 1, kjpindex
3118	DO jc = 1, nvm
3119	IF(veget_max(ic,jc) == zero)THEN
3120	! this vegetation type is not present, so no reason to do the
3121	! calculation
3122	CYCLE
3123	ENDIF
3124
3125	DO kc = 1, nlevels_tot
3126	IF (gs_column_total(ic,jc) .ne. zero) THEN
3127	gs_distribution(ic,jc,kc) = gs_store(ic,jc,kc)* lai_per_level(ic,jc,kc) / gs_column_total(ic,jc)
3128	ELSE
3129	gs_distribution(ic,jc,kc) = zero
3130	END IF
3131	Abs_light(ic,jc,kc)=Light_Abs_Tot(ic,jc,kc)*light_tran_to_level(ic,jc,kc)
3132	END DO ! kc = 1, nlevels_tot
3133	IF (jv==test_pft .AND. printlev_loc>=4) THEN
3134	Abs_light_can(:,:)=SUM(Abs_light(:,:,:),3)
3135	WRITE(numout,*) 'Absorbed light in the canopy', Abs_light_can(:,jc)
3136	ENDIF
3137	END DO ! jc = 1, nvm
3138	END DO ! ic = 1, kjpindex
3139
3140	!!$ !++++++ WRITING OUT SOME STUFF
3141	!!$ counter=counter+1
3142	!!$ WRITE(123,'(I8,20F20.10)') counter,gs_distribution(test_grid,test_pft,:)
3143	!!$ WRITE(124,'(I8,20F20.10)') counter,gstot_frac(test_grid,test_pft,:)
3144	!!$ WRITE(125,'(I8,20F20.10)') counter,gstot_store(test_grid,test_pft)
3145	!!$ WRITE(126,'(I8,20F20.10)') counter,gs_store(test_grid,test_pft,:)
3146	!!$ WRITE(127,'(I8,20F20.10)') counter,lai_per_level(test_grid,test_pft,:)
3147	!!$ WRITE(128,'(I8,20F20.10)') counter,gstot_component(test_grid,test_pft,:)
3148	!!$ WRITE(129,'(I8,20F20.10)') counter,gs_column_total(test_grid,test_pft)
3149	!!$ WRITE(130,'(I8,20F20.10)') counter,SUM(lai_per_level(test_grid,test_pft,:))
3150	!!$ IF(rveget(test_grid,test_pft) .LT. 1e19) &
3151	!!$ WRITE(131,'(I8,20ES20.10)') counter,rveget(test_grid,test_pft)
3152	!!$ IF(rveget_min(test_grid,test_pft) .LT. 1e19) &
3153	!!$ WRITE(132,'(I8,20ES20.10)') counter,rveget_min(test_grid,test_pft)
3154	!!$ WRITE(133,'(I8,20ES20.10)') counter,vbeta3(test_grid,test_pft)
3155	!!$ WRITE(134,'(I8,20ES20.10)') counter,vbeta3pot(test_grid,test_pft)
3156
3157	! Write some output for Machine learning trials
3158	CALL xios_orchidee_send_field("temp_growth",temp_growth)
3159	CALL xios_orchidee_send_field("rau",rau)
3160	CALL xios_orchidee_send_field("Ca",Ca)
3161	CALL xios_orchidee_send_field("vbeta3pot",vbeta3pot)
3162	CALL xios_orchidee_send_field("LIGHT_TRAN_TO_LEVEL",Light_Tran_Tot)
3163	CALL xios_orchidee_send_field("LIGHT_ABS_TO_LEVEL",Light_Abs_Tot)
3164
3165	IF (printlev_loc>=3) WRITE (numout,*) ' diffuco_trans_co2 done '
3166
3167	END SUBROUTINE diffuco_trans_co2
3168
3169
3170	!! ================================================================================================================================
3171	!! SUBROUTINE : diffuco_comb
3172	!!
3173	!>\BRIEF This routine combines the previous partial beta
3174	!! coefficients and calculates the total alpha and complete beta coefficients.
3175	!!
3176	!! DESCRIPTION : Those integrated coefficients are used to calculate (in enerbil.f90) the total evapotranspiration
3177	!! from the grid-cell. \n
3178	!!
3179	!! In the case that air is more humid than surface, dew deposition can occur (negative latent heat flux).
3180	!! In this instance, for temperature above zero, all of the beta coefficients are set to 0, except for
3181	!! interception (vbeta2) and bare soil (vbeta4 with zero soil resistance). The amount of water that is
3182	!! intercepted by leaves is calculated based on the value of LAI of the surface. In the case of freezing
3183	!! temperatures, water is added to the snow reservoir, and so vbeta4 and vbeta2 are set to 0, and the
3184	!! total vbeta is set to 1.\n
3185	!!
3186	!! \latexonly
3187	!! \input{diffucocomb1.tex}
3188	!! \endlatexonly
3189	!!
3190	!! The beta and alpha coefficients are initially set to 1.
3191	!! \latexonly
3192	!! \input{diffucocomb2.tex}
3193	!! \endlatexonly
3194	!!
3195	!! If snow is lower than the critical value:
3196	!! \latexonly
3197	!! \input{diffucocomb3.tex}
3198	!! \endlatexonly
3199	!! If in the presence of dew:
3200	!! \latexonly
3201	!! \input{diffucocomb4.tex}
3202	!! \endlatexonly
3203	!!
3204	!! Determine where the water goes (soil, vegetation, or snow)
3205	!! when air moisture exceeds saturation.
3206	!! \latexonly
3207	!! \input{diffucocomb5.tex}
3208	!! \endlatexonly
3209	!!
3210	!! If it is not freezing dew is put into the interception reservoir and onto the bare soil. If it is freezing,
3211	!! water is put into the snow reservoir.
3212	!! Now modify vbetas where necessary: for soil and snow
3213	!! \latexonly
3214	!! \input{diffucocomb6.tex}
3215	!! \endlatexonly
3216	!!
3217	!! and for vegetation
3218	!! \latexonly
3219	!! \input{diffucocomb7.tex}
3220	!! \endlatexonly
3221	!!
3222	!! Then compute part of dew that can be intercepted by leafs.
3223	!!
3224	!! There will be no transpiration when air moisture is too high, under any circumstance
3225	!! \latexonly
3226	!! \input{diffucocomb8.tex}
3227	!! \endlatexonly
3228	!!
3229	!! There will also be no interception loss on bare soil, under any circumstance.
3230	!! \latexonly
3231	!! \input{diffucocomb9.tex}
3232	!! \endlatexonly
3233	!!
3234	!! The flowchart details the 'decision tree' which underlies the module.
3235	!!
3236	!! RECENT CHANGE(S): None
3237	!!
3238	!! MAIN OUTPUT VARIABLE(S): vbeta1, vbeta4, humrel, vbeta2, vbeta3, vbeta
3239	!!
3240	!! REFERENCE(S) :
3241	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
3242	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
3243	!! Circulation Model. Journal of Climate, 6, pp.248-273
3244	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
3245	!! sur le cycle de l'eau, PhD Thesis, available from:
3246	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
3247	!!
3248	!! FLOWCHART :
3249	!! \latexonly
3250	!! \includegraphics[scale=0.25]{diffuco_comb_flowchart.png}
3251	!! \endlatexonly
3252	!! \n
3253	!_ ================================================================================================================================
3254
3255	SUBROUTINE diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
3256	snow, veget, lai, &
3257	tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, &
3258	evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax, vbeta2sum, vbeta3sum)
3259
3260
3261	! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006
3262
3263	!! 0. Variable and parameter declaration
3264
3265	!! 0.1 Input variables
3266
3267	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (-)
3268	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
3269	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
3270	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Nortward Lowest level wind speed (m s^{-1})
3271	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
3272	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
3273	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
3274	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
3275	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lower air temperature (K)
3276	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
3277	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (fraction)
3278	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
3279	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
3280	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil!! Total evaporating bare soil fraction
3281	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
3282
3283	!! 0.2 Output variables
3284
3285	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
3286
3287	!! 0.3 Modified variables
3288
3289	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation (-)
3290	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation (-)
3291	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
3292	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for interception loss (-)
3293	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration (-)
3294	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: vbeta2sum, vbeta3sum
3295
3296	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! limiting factor for bare soil evaporation
3297	!! when the 11-layer hydrology is used (-)
3298	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout):: evap_bare_lim_ns !! limiting factor for bare soil evaporation
3299	!! when the 11-layer hydrology is used (-)
3300	!! 0.4 Local variables
3301
3302	INTEGER(i_std) :: ji, jv
3303	REAL(r_std) :: zevtest, zsoil_moist, zrapp
3304	REAL(r_std), DIMENSION(kjpindex) :: qsatt
3305	LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow
3306	REAL(r_std) :: coeff_dew_veg
3307	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum_temp,vbeta3sum_temp
3308	REAL(r_std), DIMENSION(kjpindex) :: vegtot
3309	!_ ================================================================================================================================
3310
3311	!! \latexonly
3312	!! \input{diffucocomb1.tex}
3313	!! \endlatexonly
3314
3315	!! 0. Initialisation of vbeta2sum and vbeta3sum. They are needed in the
3316	!! calculations for the multi-layer energy budget.
3317	vbeta2sum(:) = zero
3318	vbeta3sum(:) = zero
3319
3320	!! 1. If dew is present:
3321
3322	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
3323
3324
3325	!! 1.1 Determine where the water goes
3326	!! Determine where the water goes (soil, vegetation, or snow)
3327	!! when air moisture exceeds saturation.
3328	!! \latexonly
3329	!! \input{diffucocomb5.tex}
3330	!! \endlatexonly
3331	toveg(:) = .FALSE.
3332	tosnow(:) = .FALSE.
3333	DO ji = 1, kjpindex
3334	IF ( qsatt(ji) .LT. qair(ji) ) THEN
3335	IF (temp_air(ji) .GT. tp_00) THEN
3336
3337	!! If it is not freezing dew is put into the
3338	!! interception reservoir and onto the bare soil.
3339	toveg(ji) = .TRUE.
3340	ELSE
3341
3342	!! If it is freezing water is put into the
3343	!! snow reservoir.
3344	tosnow(ji) = .TRUE.
3345	ENDIF
3346	ENDIF
3347	END DO
3348
3349	!! 1.2 Now modify vbetas where necessary.
3350
3351	!! 1.2.1 Soil and snow
3352
3353	!! \latexonly
3354	!! \input{diffucocomb6.tex}
3355	!! \endlatexonly
3356
3357	! We need to keep consistency between evap_bare_lim, evap_bare_lim_ns and vbeta4 (thus vevapnu)
3358	! or we have a water conservation issue in hydrol_split_soil
3359
3360	DO ji = 1, kjpindex
3361
3362	IF ( toveg(ji) ) THEN
3363
3364	vbeta1(ji) = zero
3365	vegtot(ji) = SUM(veget_max(ji,:))
3366
3367	IF ( (tot_bare_soil(ji) .GT. min_sechiba) .AND. (vegtot(ji).GT. min_sechiba) ) THEN
3368
3369	vbeta4(ji) = tot_bare_soil(ji)
3370
3371	evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji)
3372
3373	evap_bare_lim(ji) = vbeta4(ji)
3374	! consistent with evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)soiltile(ji,:)vegtot(ji))
3375	! as SUM(soiltile(ji,:)) = 1
3376
3377	ELSE
3378	vbeta4(ji) = zero
3379	evap_bare_lim_ns(ji,:) = zero
3380	evap_bare_lim(ji) = zero
3381	ENDIF
3382	ENDIF
3383
3384	IF ( tosnow(ji) ) THEN
3385	vbeta1(ji) = un
3386	vbeta4(ji) = zero
3387	evap_bare_lim_ns(ji,:) = zero
3388	evap_bare_lim(ji) = zero
3389	ENDIF
3390
3391	ENDDO
3392
3393	!! 1.2.2 Vegetation and interception loss
3394	!! \latexonly
3395	!! \input{diffucocomb7.tex}
3396	!! \endlatexonly
3397	DO jv = 1, nvm
3398
3399	DO ji = 1, kjpindex
3400	IF ( toveg(ji) ) THEN
3401	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
3402
3403	! Compute part of dew that can be intercepted by leafs.
3404	IF ( lai(ji,jv) .GT. min_sechiba) THEN
3405	IF (lai(ji,jv) .GT. 1.5) THEN
3406	coeff_dew_veg= &
3407	& dew_veg_poly_coeff(6)lai(ji,jv)*5 &
3408	& - dew_veg_poly_coeff(5)lai(ji,jv)*4 &
3409	& + dew_veg_poly_coeff(4)lai(ji,jv)*3 &
3410	& - dew_veg_poly_coeff(3)lai(ji,jv)*2 &
3411	& + dew_veg_poly_coeff(2)*lai(ji,jv) &
3412	& + dew_veg_poly_coeff(1)
3413	ELSE
3414	coeff_dew_veg=un
3415	ENDIF
3416	ELSE
3417	coeff_dew_veg=zero
3418	ENDIF
3419	IF (jv .EQ. 1) THEN
3420	! This line may not work with CWRR when frac_bare is distributed
3421	! among three soiltiles. Fortunately, qsintmax(ji,1)=0 (LAI=0 in
3422	! PFT1) so we never pass here
3423	vbeta2(ji,jv) = coeff_dew_veg*tot_bare_soil(ji)
3424	ELSE
3425	vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv)
3426	ENDIF
3427	ELSE
3428	vbeta2(ji,jv) = zero ! if qsintmax=0, vbeta2=0
3429	ENDIF
3430	ENDIF
3431	IF ( tosnow(ji) ) vbeta2(ji,jv) = zero
3432
3433	ENDDO
3434
3435	ENDDO
3436
3437	!! 1.2.3 Vegetation and transpiration
3438	!! There will be no transpiration when air moisture is too high, under any circumstance
3439	!! \latexonly
3440	!! \input{diffucocomb8.tex}
3441	!! \endlatexonly
3442	DO jv = 1, nvm
3443	DO ji = 1, kjpindex
3444	IF ( qsatt(ji) .LT. qair(ji) ) THEN
3445	vbeta3(ji,jv) = zero
3446	humrel(ji,jv) = zero
3447	ENDIF
3448	ENDDO
3449	ENDDO
3450
3451
3452	!! 1.2.4 Overrules 1.2.2
3453	!! There will also be no interception loss on bare soil, under any circumstance.
3454	!! \latexonly
3455	!! \input{diffucocomb9.tex}
3456	!! \endlatexonly
3457	DO ji = 1, kjpindex
3458	IF ( qsatt(ji) .LT. qair(ji) ) THEN
3459	vbeta2(ji,1) = zero
3460	ENDIF
3461	ENDDO
3462
3463	!! 2. Now calculate vbeta in all cases (the equality needs to hold for enerbil to be consistent)
3464	DO ji = 1, kjpindex
3465	vbeta(ji) = vbeta4(ji) + SUM(vbeta2(ji,:)) + SUM(vbeta3(ji,:))
3466
3467	vbeta2sum(ji) = SUM(vbeta2(ji,:))
3468	vbeta3sum(ji) = SUM(vbeta3(ji,:))
3469
3470
3471	IF (vbeta(ji) .LT. min_sechiba) THEN
3472	vbeta(ji) = zero
3473	vbeta4(ji) = zero
3474	vbeta2(ji,:)= zero
3475	vbeta3(ji,:)= zero
3476	evap_bare_lim_ns(ji,:) = zero
3477	evap_bare_lim(ji) = zero
3478	END IF
3479	ENDDO
3480
3481	CALL xios_orchidee_send_field("evap_bare_lim",evap_bare_lim)
3482	CALL xios_orchidee_send_field("evap_bare_lim_ns",evap_bare_lim_ns)
3483
3484	IF (printlev>=3) WRITE (numout,*) ' diffuco_comb done '
3485
3486	END SUBROUTINE diffuco_comb
3487
3488
3489	!! ================================================================================================================================
3490	!! SUBROUTINE : diffuco_raerod
3491	!!
3492	!>\BRIEF Computes the aerodynamic resistance, for cases in which the
3493	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ and when the flag
3494	!! 'ldq_cdrag_from_gcm' is set to TRUE
3495	!!
3496	!! DESCRIPTION : Simply computes the aerodynamic resistance, for cases in which the
3497	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient
3498	!! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine
3499	!! diffuco_aero is called instead of this one.
3500	!!
3501	!! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed:
3502	!! \latexonly
3503	!! \input{diffucoaerod1.tex}
3504	!! \endlatexonly
3505	!!
3506	!! next calculate ::raero
3507	!! \latexonly
3508	!! \input{diffucoaerod2.tex}
3509	!! \endlatexonly
3510	!!
3511	!! RECENT CHANGE(S): None
3512	!!
3513	!! MAIN OUTPUT VARIABLE(S): ::raero
3514	!!
3515	!! REFERENCE(S) :
3516	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
3517	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
3518	!! Circulation Model. Journal of Climate, 6, pp.248-273
3519	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influence de l'irrigation
3520	!! sur le cycle de l'eau, PhD Thesis, available from:
3521	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
3522	!!
3523	!! FLOWCHART : None
3524	!! \n
3525	!_ ================================================================================================================================
3526
3527	SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
3528
3529	IMPLICIT NONE
3530
3531	!! 0. Variable and parameter declaration
3532
3533	!! 0.1 Input variables
3534
3535	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
3536	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind velocity (m s^{-1})
3537	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind velocity (m s^{-1})
3538	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
3539
3540	!! 0.2 Output variables
3541
3542	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance (s m^{-1})
3543
3544	!! 0.3 Modified variables
3545
3546	!! 0.4 Local variables
3547
3548	INTEGER(i_std) :: ji !! (-)
3549	REAL(r_std) :: speed !! (m s^{-1})
3550	!_ ================================================================================================================================
3551
3552	!! 1. Simple calculation of the aerodynamic resistance, for diganostic purposes.
3553
3554	DO ji=1,kjpindex
3555
3556	!! \latexonly
3557	!! \input{diffucoaerod1.tex}
3558	!! \endlatexonly
3559	speed = MAX(min_wind, wind(ji))
3560
3561	!! \latexonly
3562	!! \input{diffucoaerod2.tex}
3563	!! \endlatexonly
3564	raero(ji) = un / (q_cdrag(ji)*speed)
3565
3566	ENDDO
3567
3568	END SUBROUTINE diffuco_raerod
3569
3570
3571	!! ================================================================================================================================
3572	!! SUBROUTINE : isotope_c13
3573	!!
3574	!>\BRIEF Calculate the fractionation of C13 during photosynthesis\n
3575	!!
3576	!! DESCRIPTION :
3577	!!
3578	!! RECENT CHANGE(S): None
3579	!!
3580	!! MAIN OUTPUT VARIABLE(S): :: delta_c13_assim
3581	!!
3582	!! REFERENCE(S) :
3583	!! - Farquhar, G. D., Ehleringer, J. R., & Hubick, K. T. (1989). CARBON ISOTOPE DISCRIMINATION AND PHOTOSYNTHESIS.
3584	!! Annual Review of Plant Physiology and Plant Molecular Biology, 40, 503-537. doi: 10.1146/annurev.arplant.40.1.50
3585	!!
3586	!! FLOWCHART : None
3587	!_ ================================================================================================================================
3588
3589	SUBROUTINE isotope_c13(kjpindex, index_assi, leaf_ci, ccanopy, lai_per_level, &
3590	assimtot, delta_c13_assim, leaf_ci_out, nvm, &
3591	nlevels_tot, assimi_lev)
3592
3593	IMPLICIT NONE
3594
3595	!! 0. Variables and parameter declaration
3596
3597	!! 0.1 Input variables
3598	INTEGER(i_std), INTENT (in) :: kjpindex !! number of land pixels
3599	INTEGER(i_std), INTENT (in) :: nvm !! number of pfts
3600	INTEGER(i_std), INTENT (in) :: nlevels_tot !! number of lai layers in canopy
3601	INTEGER(i_std), DIMENSION(:), INTENT (in) :: index_assi !! indices (unitless)
3602	REAL(r_std), DIMENSION (:,:,:), INTENT (in) :: leaf_ci !! intercellular CO2 concentration (ppm)
3603	REAL(r_std),DIMENSION (:), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
3604	REAL(r_std), DIMENSION(:,:,:), INTENT (in) :: lai_per_level !! This is the LAI per vertical level
3605	!! @tex $(m^{2} m^{-2})$ @endtex
3606	!+++CHECK+++
3607	! Looks like assimi is a local variable in trans_co2 and that
3608	! its values are overwritten for every new layer. If this is
3609	! indeed the case, it should not be used outside trans_co2
3610	! REAL(r_std), DIMENSION(:,:), INTENT (in) :: assimi !! assimilation (at a specific LAI level)
3611	! !! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
3612	!+++++++++++
3613	REAL(r_std), DIMENSION(:,:), INTENT (in) :: assimtot !! total assimilation
3614	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
3615	REAL(r_std), DIMENSION(:,:,:), INTENT (in) :: assimi_lev !! assimilation (at all LAI levels)
3616	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
3617
3618	!! 0.2 Output variables
3619	REAL(r_std),DIMENSION (:,:), INTENT (out) :: delta_c13_assim !! C13 concentration in delta notation
3620	!! @tex $ permille $ @endtex (per thousand)
3621	REAL(r_std),DIMENSION (:,:), INTENT (out) :: leaf_ci_out !! Ci @tex $ ??unit?? $ @endtex
3622
3623	!! 0.3 Modified variables
3624
3625	!! 0.4 Local variables
3626	INTEGER(i_std) :: inia,iainia !! counter/indices (unitless)
3627	INTEGER(i_std) :: ji,jv,jl !! Indices (-)
3628
3629
3630	!_ ================================================================================================================================
3631
3632	!! Initialize
3633	delta_c13_assim(:,:) = zero
3634	leaf_ci_out(:,:) = zero
3635
3636	!! Calculate delta_C13
3637	DO ji=1,kjpindex
3638
3639	! Skip the whole pixel of none of the PFTs did photosynthesis
3640	IF (SUM(SUM(assimi_lev(ji,:,:),2),1) .LE. zero) THEN
3641	! There was no photsynthesis at this pixel, no reason to do the
3642	! calculation
3643	CYCLE
3644	ENDIF
3645
3646	DO jv = 1,nvm
3647
3648	! Skip the PFT is it didn't do photosynthesis
3649	IF (SUM(assimi_lev(ji,jv,:)) .LE. zero) THEN
3650	! There was no photsynthesis at this pixel, no reason to do the
3651	! calculation
3652	CYCLE
3653	ENDIF
3654
3655	! There was photosynthesis at this PFT, calculate the C13 fractionation
3656	DO jl = 1,nlevels_tot
3657
3658	IF (lai_per_level(ji,jv,jl) .GT. zero) THEN
3659
3660	IF (printlev_loc>=4) THEN
3661	WRITE(numout,*) 'leaf_ci, ',leaf_ci(ji,jv,jl)
3662	WRITE(numout,*) 'assimi_lev, ',assimi_lev(ji,jv,jl)
3663	WRITE(numout,*) 'assimtot - isotope - before calc, ',assimtot(ji,jv)
3664	WRITE(numout,*) 'lai_per_level, ',lai_per_level(ji,jv,jl)
3665	WRITE(numout,*) 'day, ',day_end
3666	WRITE(numout,*) 'pft, ',jv
3667	ENDIF
3668
3669	! Simple version of Farquhar model
3670	! delta_c13_assim(:,jv) = (-8 - c13_a - (c13_b - c13_a) * cimean(:,jv) / ccanopy(:))
3671	delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv) + (- c13_a - (c13_b - c13_a)* &
3672	leaf_ci(ji,jv,jl)/ccanopy(ji)) * &
3673	assimi_lev(ji,jv,jl) * lai_per_level(ji,jv,jl)
3674	leaf_ci_out(ji,jv) = leaf_ci_out(ji,jv) + &
3675	leaf_ci(ji,jv,jl) * &
3676	assimi_lev(ji,jv,jl) * lai_per_level(ji,jv,jl)
3677
3678	ENDIF
3679
3680	! Simple version of Farquhar model with photorespiratioin, gi and gamma*
3681	! This code copied from the work by Thomas Eglin and Thomas Launois where it
3682	! was already commented out. When moving it to ORCHIDEE-CAN is has not been
3683	! tested with ORCHIDEE-CAN. It is not clear why it was commented out.
3684	! was commented out
3685	!!$ IF (assimi(ji,jv) .GT. 0) THEN
3686	!!$ delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv) + (&
3687	!!$ - c13_ab * (ccanopy(ji)-leaf_cs(ji,jv,jl))/(ccanopy(ji)) &
3688	!!$ - c13_a * (leaf_cs(ji,jv,jl)-leaf_ci(ji,jv,jl))/(ccanopy(ji)) &
3689	!!$ - (C13_es + c13_al) * (leaf_ci(ji,jv,jl)-leaf_cc(ji,jv,jl))/(ccanopy(ji)) &
3690	!!$ - (c13_b+2.0) * leaf_cc(ji,jv,jl)/(ccanopy(ji)) &
3691	!!$ - (C13_e * (Rd(ji)/(Rd(ji) + assimi(ji,jv))) + C13_f * CP(ji))/(ccanopy(ji))&
3692	!!$ ) * assimi(ji,jv) * (laitab(jl+1)-laitab(jl))
3693	!!$ ENDIF
3694
3695
3696	!!$ (Rd(ji) + assimi(ji,jv)), Rd(ji),assimi(ji,jv)
3697
3698	! This code copied from the work by Thomas Eglin and Thomas Launois where it
3699	! was already commented out. When moving it to ORCHIDEE-CAN is has not been
3700	! tested with ORCHIDEE-CAN. It is not clear why it was commented out.
3701	! was commented out
3702	! Modification M; Bordigoni 07/2008 : New calculation of c3-plants fractionation
3703	! according to Lloyd and Farquhar_Oecologia_1994
3704	!!$ delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv) +(- (&
3705	!!$ 4.4*(1-(leaf_ci(ji,jv,jl) / ccanopy(ji) )+0.025)+ &
3706	!!$ 0.075*(1.1+0.7)+ &
3707	!!$ 27.5*(( leaf_ci(ji,jv,jl) / ccanopy(ji))-0.1)- &
3708	!!$ (110.1+8(1.541.05((temp_air(ji)-273.16)+2.5)))/ ccanopy(ji)) )* &
3709	!!$ assimi(ji,jv) * (laitab(jl+1)-laitab(jl))
3710
3711	ENDDO ! photosynthesis layers
3712
3713	IF (assimtot(ji,jv) .GT. threshold_c13_assim) THEN
3714
3715	delta_c13_assim(ji,jv) = delta_c13_assim(ji,jv)/assimtot(ji,jv)
3716	leaf_ci_out(ji,jv) = leaf_ci_out(ji,jv)/assimtot(ji,jv)
3717
3718	ELSE
3719
3720	delta_c13_assim(ji,jv) = zero
3721	leaf_ci_out(ji,jv) = zero
3722
3723	ENDIF
3724
3725	ENDDO ! # PFTs
3726
3727	ENDDO ! #kjpindex
3728
3729	END SUBROUTINE isotope_c13
3730
3731
3732	END MODULE diffuco

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

Download in other formats: