Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/diffuco.f90 @ 7485

Last change on this file since 7485 was 7326, checked in by josefine.ghattas, 3 years ago
Corrected bug on carbon balance closure. See ticket #785 Integration in branch 2_2 done by P. Cadule
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 134.5 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
61	USE ioipsl_para
62	USE xios_orchidee
63	USE chemistry, ONLY : chemistry_initialize, chemistry_bvoc, chemistry_clear
64	IMPLICIT NONE
65
66	! public routines :
67	PRIVATE
68	PUBLIC :: diffuco_main, diffuco_initialize, diffuco_finalize, diffuco_clear
69
70	INTERFACE Arrhenius_modified
71	MODULE PROCEDURE Arrhenius_modified_0d, Arrhenius_modified_1d
72	END INTERFACE
73
74	!
75	! variables used inside diffuco module : declaration and initialisation
76	!
77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind module (m s^{-1})
78	!$OMP THREADPRIVATE(wind)
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	!! 1. Define flag ldq_cdrag_from_gcm. This flag determines if the cdrag should be taken from the GCM or be calculated.
127	!! The default value is true if the q_cdrag variables was already initialized. This is the case when coupled to the LMDZ.
128
129	!Config Key = CDRAG_FROM_GCM
130	!Config Desc = Keep cdrag coefficient from gcm.
131	!Config If = OK_SECHIBA
132	!Config Def = y
133	!Config Help = Set to .TRUE. if you want q_cdrag coming from GCM (if q_cdrag on initialization is non zero).
134	!Config Keep cdrag coefficient from gcm for latent and sensible heat fluxes.
135	!Config Units = [FLAG]
136	IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN
137	ldq_cdrag_from_gcm = .FALSE.
138	ELSE
139	ldq_cdrag_from_gcm = .TRUE.
140	ENDIF
141	CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm)
142	IF (printlev>=2) WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm
143
144	!! 2. Allocate module variables
145	ALLOCATE (wind(kjpindex),stat=ier)
146	IF (ier /= 0) CALL ipslerr_p(3,'diffuco_initialize','Problem in allocate of variable wind','','')
147
148	!! 3. Read variables from restart file
149	IF (printlev>=3) WRITE (numout,*) 'Read DIFFUCO variables from restart file'
150
151	CALL ioconf_setatt_p('UNITS', 's/m')
152	CALL ioconf_setatt_p('LONG_NAME','Structural resistance')
153	CALL restget_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g)
154	IF ( ALL(rstruct(:,:) == val_exp) ) THEN
155	DO jv = 1, nvm
156	rstruct(:,jv) = rstruct_const(jv)
157	ENDDO
158	ENDIF
159
160	!! 4. Initialize chemistry module
161	IF (printlev>=3) WRITE(numout,*) "ok_bvoc:",ok_bvoc
162	IF ( ok_bvoc ) CALL chemistry_initialize(kjpindex, lalo, neighbours, resolution)
163
164	END SUBROUTINE diffuco_initialize
165
166
167
168	!! ================================================================================================================================
169	!! SUBROUTINE : diffuco_main
170	!!
171	!>\BRIEF The root subroutine for the module, which calls all other required
172	!! subroutines.
173	!!
174	!! DESCRIPTION :
175
176	!! This is the main subroutine for the module.
177	!! First it calculates the surface drag coefficient (via a call to diffuco_aero), using available parameters to determine
178	!! the stability of air in the surface layer by calculating the Richardson Nubmber. If a value for the
179	!! surface drag coefficient is passed down from the atmospheric model and and if the flag 'ldq_cdrag_from_gcm'
180	!! is set to TRUE, then the subroutine diffuco_aerod is called instead. This calculates the aerodynamic coefficient. \n
181	!!
182	!! Following this, an estimation of the saturated humidity at the surface is made (via a call
183	!! to qsatcalc in the module qsat_moisture). Following this the beta coefficients for sublimation (via
184	!! diffuco_snow), interception (diffuco_inter), bare soil (diffuco_bare), and transpiration (via
185	!! diffuco_trans_co2) are calculated in sequence. Finally
186	!! the alpha and beta coefficients are combined (diffuco_comb). \n
187	!!
188	!! The surface drag coefficient is calculated for use within the module enerbil. It is required to to
189	!! calculate the aerodynamic coefficient for every flux. \n
190	!!
191	!! The various beta coefficients are used within the module enerbil for modifying the rate of evaporation,
192	!! as appropriate for the surface. As explained in Chapter 2 of Guimberteau (2010), that module (enerbil)
193	!! calculates the rate of evaporation essentially according to the expression $E = /beta E_{pot}$, where
194	!! E is the total evaporation and $E_{pot}$ the evaporation potential. If $\beta = 1$, there would be
195	!! essentially no resistance to evaporation, whereas should $\beta = 0$, there would be no evaporation and
196	!! the surface layer would be subject to some very stong hydrological stress. \n
197	!!
198	!! The following processes are calculated:
199	!! - call diffuco_aero for aerodynamic transfer coeficient
200	!! - call diffuco_snow for partial beta coefficient: sublimation
201	!! - call diffuco_inter for partial beta coefficient: interception for each type of vegetation
202	!! - call diffuco_bare for partial beta coefficient: bare soil
203	!! - call diffuco_trans_co2 for partial beta coefficient: transpiration for each type of vegetation, using Farqhar's formula
204	!! - call diffuco_comb for alpha and beta coefficient
205	!! - call chemistry_bvoc for alpha and beta coefficients for biogenic emissions
206	!!
207	!! RECENT CHANGE(S): None
208	!!
209	!! MAIN OUTPUT VARIABLE(S): humrel, q_cdrag, vbeta, vbeta1, vbeta4,
210	!! vbeta2, vbeta3, rveget, cimean
211	!!
212	!! REFERENCE(S) :
213	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
214	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
215	!! Circulation Model. Journal of Climate, 6, pp.248-273.
216	!! - de Rosnay, P, 1999. ReprÃ©sentation des interactions sol-plante-atmosphÃšre dans le modÃšle de circulation gÃ©nÃ©rale
217	!! du LMD, 1999. PhD Thesis, UniversitÃ© Paris 6, available (25/01/12):
218	!! http://www.ecmwf.int/staff/patricia_de_rosnay/publications.html#8
219	!! - Ducharne, A, 1997. Le cycle de l'eau: modÃ©lisation de l'hydrologie continentale, Ã©tude de ses interactions avec
220	!! le climat, PhD Thesis, UniversitÃ© Paris 6
221	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
222	!! sur le cycle de l'eau, PhD Thesis, available (25/01/12):
223	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
224	!! - LathiÃšre, J, 2005. Evolution des Ã©missions de composÃ©s organiques et azotÃ©s par la biosphÃšre continentale dans le
225	!! modÃšle LMDz-INCA-ORCHIDEE, UniversitÃ© Paris 6
226	!!
227	!! FLOWCHART :
228	!! \latexonly
229	!! \includegraphics[scale=0.5]{diffuco_main_flowchart.png}
230	!! \endlatexonly
231	!! \n
232	!_ ================================================================================================================================
233
234	SUBROUTINE diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
235	& zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, q_cdrag, qsurf, qair, pb, &
236	& evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
237	& swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
238	& vbeta , vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, &
239	& lalo, neighbours, resolution, ptnlev1, precip_rain, frac_age, tot_bare_soil, frac_snow_veg, frac_snow_nobio, &
240	& hist_id, hist2_id)
241
242	!! 0. Variable and parameter declaration
243
244	!! 0.1 Input variables
245
246	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
247	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
248	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-)
249	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (-)
250	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
251	INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai !! Indeces of the points on the 3D map
252	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-)
253	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
254	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
255	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
256	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m)
257	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m)
258	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness (m)
259	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
260	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
261	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
262	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
263	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface air specific humidity (kg kg^{-1})
264	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level air specific humidity (kg kg^{-1})
265	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
266	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! Fraction of floodplains
267	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! Reservoir in floodplains (estimation to avoid over-evaporation)
268	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure (hPa)
269	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! Limit to the bare soil evaporation when the
270	!! 11-layer hydrology is used (-)
271	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! Limit to the bare soil evaporation when the
272	!! 11-layer hydrology is used (-)
273	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1})
274	!! NdN This variable does not seem to be used at
275	!! all in diffuco
276	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation
277	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... (-)
278	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... (kg m^{-2})
279	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... (-)
280	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
281	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux (W m^{-2})
282	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: coszang !! Cosine of the solar zenith angle (unitless)
283	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
284	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (-)
285	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
286	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
287	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
288	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
289	!! (kg m^{-2})
290	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! min+max+opt temps, vcmax, vjmax
291	!! for photosynthesis (K ??)
292	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
293	INTEGER(i_std),DIMENSION (kjpindex,NbNeighb),INTENT (in):: neighbours !! Vector of neighbours for each
294	!! grid point (1=N, 2=E, 3=S, 4=W)
295	REAL(r_std),DIMENSION (kjpindex,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
296	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ptnlev1 !! 1st level of soil temperature (Kelvin)
297	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation expressed in mm/tstep
298	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (in) :: frac_age !! Age efficiency for isoprene emissions (from STOMATE)
299	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
300	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
301	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in):: frac_snow_nobio !! Snow cover fraction on non-vegeted area
302
303	!! 0.2 Output variables
304
305	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
306	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (-)
307	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
308	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains
309	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! Mean stomatal conductance to CO2 (mol m-2 s-1)
310	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
311	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration (-)
312	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for potential transpiration
313	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Stomatal resistance for the whole canopy (s m^{-1})
314	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! Structural resistance for the vegetation
315	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! Mean leaf Ci (ppm)
316	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: gpp !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area)
317
318	!! 0.3 Modified variables
319
320	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
321	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag coefficient (-)
322
323	!! 0.4 Local variables
324
325	REAL(r_std),DIMENSION (kjpindex,nvm) :: vbeta23 !! Beta for fraction of wetted foliage that will
326	!! transpire once intercepted water has evaporated (-)
327	REAL(r_std),DIMENSION (kjpindex) :: raero !! Aerodynamic resistance (s m^{-1})
328	INTEGER(i_std) :: ilai
329	CHARACTER(LEN=4) :: laistring
330	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
331	REAL(r_std),DIMENSION(kjpindex) :: qsatt !! Surface saturated humidity (kg kg^{-1})
332	REAL(r_std),DIMENSION (kjpindex,nvm) :: cim !! Intercellular CO2 over nlai
333
334	!_ ================================================================================================================================
335
336	wind(:) = SQRT (u(:)u(:) + v(:)v(:))
337
338	!! 1. Calculate the different coefficients
339
340	IF (.NOT.ldq_cdrag_from_gcm) THEN
341	! Case 1a)
342	CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
343	qsurf, qair, snow, q_cdrag)
344	ENDIF
345
346	! Case 1b)
347	CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
348
349	!! 2. Make an estimation of the saturated humidity at the surface
350
351	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
352
353	!! 3. Calculate the beta coefficient for sublimation
354
355	CALL diffuco_snow (kjpindex, qair, qsatt, rau, u, v, q_cdrag, &
356	snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, &
357	vbeta1)
358
359
360	CALL diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
361	& flood_frac, flood_res, vbeta5)
362
363	!! 4. Calculate the beta coefficient for interception
364
365	CALL diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
366	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
367
368
369	!! 5. Calculate the beta coefficient for transpiration
370
371	CALL diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, &
372	assim_param, ccanopy, &
373	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, &
374	rveget, rstruct, cimean, gsmean, gpp, &
375	vbeta23, hist_id, indexveg, indexlai, index, kjit, cim)
376
377	!
378	!biogenic emissions
379	!
380	IF ( ok_bvoc ) THEN
381	CALL chemistry_bvoc (kjpindex, swdown, coszang, temp_air, &
382	temp_sol, ptnlev1, precip_rain, humrel, veget_max, &
383	lai, frac_age, lalo, ccanopy, cim, wind, snow, &
384	veget, hist_id, hist2_id, kjit, index, &
385	indexlai, indexveg)
386	ENDIF
387	!
388	! combination of coefficient : alpha and beta coefficient
389	! beta coefficient for bare soil
390	!
391
392	CALL diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)
393
394	!! 6. Combine the alpha and beta coefficients
395
396	! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006
397	CALL diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, &
398	& veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, &
399	& evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax)
400
401	CALL xios_orchidee_send_field("q_cdrag",q_cdrag)
402	CALL xios_orchidee_send_field("raero",raero)
403	CALL xios_orchidee_send_field("wind",wind)
404	CALL xios_orchidee_send_field("qsatt",qsatt)
405	CALL xios_orchidee_send_field("coszang",coszang)
406	CALL xios_orchidee_send_field('cim', cim)
407
408	IF ( .NOT. almaoutput ) THEN
409	CALL histwrite_p(hist_id, 'raero', kjit, raero, kjpindex, index)
410	CALL histwrite_p(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
411	CALL histwrite_p(hist_id, 'Wind', kjit, wind, kjpindex, index)
412	CALL histwrite_p(hist_id, 'qsatt', kjit, qsatt, kjpindex, index)
413	CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)
414
415	IF ( hist2_id > 0 ) THEN
416	CALL histwrite_p(hist2_id, 'raero', kjit, raero, kjpindex, index)
417	CALL histwrite_p(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
418	CALL histwrite_p(hist2_id, 'Wind', kjit, wind, kjpindex, index)
419	CALL histwrite_p(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index)
420	ENDIF
421	ELSE
422	CALL histwrite_p(hist_id, 'cim', kjit, cim, kjpindex*nvm, indexveg)
423	ENDIF
424
425	IF (printlev>=3) WRITE (numout,*) ' diffuco_main done '
426
427	END SUBROUTINE diffuco_main
428
429	!! =============================================================================================================================
430	!! SUBROUTINE: diffuco_finalize
431	!!
432	!>\BRIEF Write to restart file
433	!!
434	!! DESCRIPTION: This subroutine writes the module variables and variables calculated in diffuco
435	!! to restart file
436	!!
437	!! RECENT CHANGE(S): None
438	!! REFERENCE(S): None
439	!! FLOWCHART: None
440	!! \n
441	!_ ==============================================================================================================================
442	SUBROUTINE diffuco_finalize (kjit, kjpindex, rest_id, rstruct )
443
444	!! 0. Variable and parameter declaration
445	!! 0.1 Input variables
446	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
447	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
448	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-)
449	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! Structural resistance for the vegetation
450
451	!! 0.4 Local variables
452	INTEGER :: ilai
453	CHARACTER(LEN=4) :: laistring
454	CHARACTER(LEN=80) :: var_name
455
456	!_ ================================================================================================================================
457
458	!! 1. Prepare the restart file for the next simulation
459	IF (printlev>=3) WRITE (numout,*) 'Complete restart file with DIFFUCO variables '
460
461	CALL restput_p (rest_id, 'rstruct', nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g)
462
463	END SUBROUTINE diffuco_finalize
464
465
466	!! ================================================================================================================================
467	!! SUBROUTINE : diffuco_clear
468	!!
469	!>\BRIEF Housekeeping module to deallocate the variables
470	!! rstruct and raero
471	!!
472	!! DESCRIPTION : Housekeeping module to deallocate the variables
473	!! rstruct and raero
474	!!
475	!! RECENT CHANGE(S) : None
476	!!
477	!! MAIN OUTPUT VARIABLE(S) : None
478	!!
479	!! REFERENCE(S) : None
480	!!
481	!! FLOWCHART : None
482	!! \n
483	!_ ================================================================================================================================
484
485	SUBROUTINE diffuco_clear()
486
487	! Deallocate and reset variables in chemistry module
488	CALL chemistry_clear
489
490	END SUBROUTINE diffuco_clear
491
492
493	!! ================================================================================================================================
494	!! SUBROUTINE : diffuco_aero
495	!!
496	!>\BRIEF This module first calculates the surface drag
497	!! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled
498	!! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE
499	!!
500	!! DESCRIPTION : Computes the surface drag coefficient, for cases
501	!! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the
502	!! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric
503	!! stability in the surface layer. The formulation used to find this surface drag coefficient is
504	!! dependent on the stability determined. \n
505	!!
506	!! Designation of wind speed
507	!! \latexonly
508	!! \input{diffucoaero1.tex}
509	!! \endlatexonly
510	!!
511	!! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson
512	!! eqn.4.47, 2005). (required for calculation of the Richardson Number)
513	!! \latexonly
514	!! \input{diffucoaero2.tex}
515	!! \endlatexonly
516	!!
517	!! \latexonly
518	!! \input{diffucoaero3.tex}
519	!! \endlatexonly
520	!!
521	!! Calculation of the virtual air temperature at the surface (required for calculation
522	!! of the Richardson Number)
523	!! \latexonly
524	!! \input{diffucoaero4.tex}
525	!! \endlatexonly
526	!!
527	!! Calculation of the virtual surface temperature (required for calculation of th
528	!! Richardson Number)
529	!! \latexonly
530	!! \input{diffucoaero5.tex}
531	!! \endlatexonly
532	!!
533	!! Calculation of the squared wind shear (required for calculation of the Richardson
534	!! Number)
535	!! \latexonly
536	!! \input{diffucoaero6.tex}
537	!! \endlatexonly
538	!!
539	!! Calculation of the Richardson Number. The Richardson Number is defined as the ratio
540	!! of potential to kinetic energy, or, in the context of atmospheric science, of the
541	!! generation of energy by wind shear against consumption
542	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
543	!! becomes turbulent and vise versa). It is approximated using the expression below:
544	!! \latexonly
545	!! \input{diffucoaero7.tex}
546	!! \endlatexonly
547	!!
548	!! The Richardson Number hence calculated is subject to a minimum value:
549	!! \latexonly
550	!! \input{diffucoaero8.tex}
551	!! \endlatexonly
552	!!
553	!! Computing the drag coefficient. We add the add the height of the vegetation to the
554	!! level height to take into account that the level 'seen' by the vegetation is actually
555	!! the top of the vegetation. Then we we can subtract the displacement height.
556	!! \latexonly
557	!! \input{diffucoaero9.tex}
558	!! \endlatexonly
559	!!
560	!! For the stable case (i.e $R_i$ $\geq$ 0)
561	!! \latexonly
562	!! \input{diffucoaero10.tex}
563	!! \endlatexonly
564	!!
565	!! \latexonly
566	!! \input{diffucoaero11.tex}
567	!! \endlatexonly
568	!!
569	!! For the unstable case (i.e. $R_i$ < 0)
570	!! \latexonly
571	!! \input{diffucoaero12.tex}
572	!! \endlatexonly
573	!!
574	!! \latexonly
575	!! \input{diffucoaero13.tex}
576	!! \endlatexonly
577	!!
578	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
579	!! To prevent this, a minimum limit to the drag coefficient is defined as:
580	!!
581	!! \latexonly
582	!! \input{diffucoaero14.tex}
583	!! \endlatexonly
584	!!
585	!! RECENT CHANGE(S): None
586	!!
587	!! MAIN OUTPUT VARIABLE(S): q_cdrag
588	!!
589	!! REFERENCE(S) :
590	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
591	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
592	!! Circulation Model. Journal of Climate, 6, pp.248-273
593	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
594	!! sur le cycle de l'eau, PhD Thesis, available from:
595	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
596	!! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge
597	!! University Press, ISBN 0-521-54865-9
598	!!
599	!! FLOWCHART :
600	!! \latexonly
601	!! \includegraphics[scale=0.5]{diffuco_aero_flowchart.png}
602	!! \endlatexonly
603	!! \n
604	!_ ================================================================================================================================
605
606	SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0h, z0m, roughheight, temp_sol, temp_air, &
607	qsurf, qair, snow, q_cdrag)
608
609	!! 0. Variable and parameter declaration
610
611	!! 0.1 Input variables
612
613	INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size
614	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
615	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
616	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first atmospheric layer (m)
617	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0h !! Surface roughness Length for heat (m)
618	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0m !! Surface roughness Length for momentum (m)
619	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height (m)
620	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Ground temperature (K)
621	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
622	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific air humidity (kg kg^{-1})
623	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
624	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
625
626	!! 0.2 Output variables
627
628	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: q_cdrag !! Surface drag coefficient (-)
629
630	!! 0.3 Modified variables
631
632	!! 0.4 Local variables
633
634	INTEGER(i_std) :: ji, jv
635	REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2
636	REAL(r_std) :: zri, cd_neut, zscf, cd_tmp
637	!_ ================================================================================================================================
638
639	!! 1. Initialisation
640
641	! test if we have to work with q_cdrag or to calcul it
642	DO ji=1,kjpindex
643
644	!! 1a).1 Designation of wind speed
645	!! \latexonly
646	!! \input{diffucoaero1.tex}
647	!! \endlatexonly
648	speed = wind(ji)
649
650	!! 1a).2 Calculation of geopotentiel
651	!! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required
652	!! for calculation of the Richardson Number)
653	!! \latexonly
654	!! \input{diffucoaero2.tex}
655	!! \endlatexonly
656	zg = zlev(ji) * cte_grav
657
658	!! \latexonly
659	!! \input{diffucoaero3.tex}
660	!! \endlatexonly
661	zdphi = zg/cp_air
662
663	!! 1a).3 Calculation of the virtual air temperature at the surface
664	!! required for calculation of the Richardson Number
665	!! \latexonly
666	!! \input{diffucoaero4.tex}
667	!! \endlatexonly
668	ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji))
669
670	!! 1a).4 Calculation of the virtual surface temperature
671	!! required for calculation of the Richardson Number
672	!! \latexonly
673	!! \input{diffucoaero5.tex}
674	!! \endlatexonly
675	ztvs = temp_sol(ji) * (un + retv * qsurf(ji))
676
677	!! 1a).5 Calculation of the squared wind shear
678	!! required for calculation of the Richardson Number
679	!! \latexonly
680	!! \input{diffucoaero6.tex}
681	!! \endlatexonly
682	zdu2 = MAX(cepdu2,speed**2)
683
684	!! 1a).6 Calculation of the Richardson Number
685	!! The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the
686	!! context of atmospheric science, of the generation of energy by wind shear against consumption
687	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
688	!! becomes turbulent and vise versa).\n
689	!! It is approximated using the expression below:
690	!! \latexonly
691	!! \input{diffucoaero7.tex}
692	!! \endlatexonly
693	zri = zg * (ztvd - ztvs) / (zdu2 * ztvd)
694
695	!! The Richardson Number hence calculated is subject to a minimum value:
696	!! \latexonly
697	!! \input{diffucoaero8.tex}
698	!! \endlatexonly
699	zri = MAX(MIN(zri,5.),-5.)
700
701	!! 1a).7 Computing the drag coefficient
702	!! We add the add the height of the vegetation to the level height to take into account
703	!! that the level 'seen' by the vegetation is actually the top of the vegetation. Then we
704	!! we can subtract the displacement height.
705	!! \latexonly
706	!! \input{diffucoaero9.tex}
707	!! \endlatexonly
708
709	!! 7.0 Snow smoothering
710	!! Snow induces low levels of turbulence.
711	!! Sensible heat fluxes can therefore be reduced of ~1/3. Pomeroy et al., 1998
712	cd_neut = ct_karman ** 2. / ( LOG( (zlev(ji) + roughheight(ji)) / z0m(ji) ) * LOG( (zlev(ji) + roughheight(ji)) / z0h(ji) ) )
713
714	!! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0)
715	IF (zri .GE. zero) THEN
716
717	!! \latexonly
718	!! \input{diffucoaero10.tex}
719	!! \endlatexonly
720	zscf = SQRT(un + cd * ABS(zri))
721
722	!! \latexonly
723	!! \input{diffucoaero11.tex}
724	!! \endlatexonly
725	cd_tmp=cd_neut/(un + trois * cb * zri * zscf)
726	ELSE
727
728	!! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0)
729	!! \latexonly
730	!! \input{diffucoaero12.tex}
731	!! \endlatexonly
732	zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * &
733	& ((zlev(ji) + roughheight(ji)) / z0m(ji))))
734
735	!! \latexonly
736	!! \input{diffucoaero13.tex}
737	!! \endlatexonly
738	cd_tmp=cd_neut * (un - trois * cb * zri * zscf)
739	ENDIF
740
741	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
742	!! To prevent this, a minimum limit to the drag coefficient is defined as:
743
744	!! \latexonly
745	!! \input{diffucoaero14.tex}
746	!! \endlatexonly
747	!!
748	q_cdrag(ji) = MAX(cd_tmp, min_qc/MAX(speed,min_wind))
749
750	! In some situations it might be useful to give an upper limit on the cdrag as well.
751	! The line here should then be uncommented.
752	!q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind))
753
754	END DO
755
756	IF (printlev>=3) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done '
757
758	END SUBROUTINE diffuco_aero
759
760
761	!! ================================================================================================================================
762	!! SUBROUTINE : diffuco_snow
763	!!
764	!>\BRIEF This subroutine computes the beta coefficient for snow sublimation.
765	!!
766	!! DESCRIPTION : This routine computes beta coefficient for snow sublimation, which
767	!! integrates the snow on both vegetation and other surface types (e.g. ice, lakes,
768	!! cities etc.) \n
769	!!
770	!! A critical depth of snow (snowcri) is defined to calculate the fraction of each grid-cell
771	!! that is covered with snow (snow/snowcri) while the remaining part is snow-free.
772	!! We also carry out a first calculation of sublimation (subtest) to lower down the beta
773	!! coefficient if necessary (if subtest > snow). This is a predictor-corrector test.
774	!!
775	!! RECENT CHANGE(S): None
776	!!
777	!! MAIN OUTPUT VARIABLE(S): ::vbeta1
778	!!
779	!! REFERENCE(S) :
780	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
781	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
782	!! Circulation Model. Journal of Climate, 6, pp. 248-273
783	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
784	!! sur le cycle de l'eau, PhD Thesis, available from:
785	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
786	!!
787	!! FLOWCHART : None
788	!! \n
789	!_ ================================================================================================================================
790
791	SUBROUTINE diffuco_snow (kjpindex, qair, qsatt, rau, u, v,q_cdrag, &
792	& snow, frac_nobio, totfrac_nobio, snow_nobio, frac_snow_veg, frac_snow_nobio, &
793	vbeta1)
794
795	!! 0. Variable and parameter declaration
796
797	!! 0.1 Input variables
798
799	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
800	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
801	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
802	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
803	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
804	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
805	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
806	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg m^{-2})
807	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. (-)
808	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice, lakes, cities etc. (-)
809	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice, lakes, cities etc. (-)
810	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: frac_snow_veg !! Snow cover fraction on vegeted area
811	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_snow_nobio!! Snow cover fraction on non-vegeted area
812
813	!! 0.2 Output variables
814
815	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (dimensionless ratio)
816
817	!! 0.3 Modified variables
818
819	!! 0.4 Local variables
820
821	REAL(r_std) :: subtest !! Sublimation for test (kg m^{-2})
822	REAL(r_std) :: zrapp !! Modified factor (ratio)
823	REAL(r_std) :: speed !! Wind speed (m s^{-1})
824	REAL(r_std) :: vbeta1_add !! Beta for sublimation (ratio)
825	INTEGER(i_std) :: ji, jv !! Indices (-)
826	!_ ================================================================================================================================
827
828	!! 1. Calculate beta coefficient for snow sublimation on the vegetation\n
829
830	DO ji=1,kjpindex ! Loop over # pixels - domain size
831
832	! Fraction of mesh that can sublimate snow
833	vbeta1(ji) = (un - totfrac_nobio(ji)) * frac_snow_veg(ji)
834
835	! Limitation of sublimation in case of snow amounts smaller than the atmospheric demand.
836	speed = MAX(min_wind, wind(ji))
837
838	subtest = dt_sechiba * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * &
839	& ( qsatt(ji) - qair(ji) )
840
841	IF ( subtest .GT. min_sechiba ) THEN
842	zrapp = snow(ji) / subtest
843	IF ( zrapp .LT. un ) THEN
844	vbeta1(ji) = vbeta1(ji) * zrapp
845	ENDIF
846	ENDIF
847
848	END DO ! Loop over # pixels - domain size
849
850	!! 2. Add the beta coefficients calculated from other surfaces types (snow on ice,lakes, cities...)
851
852	DO jv = 1, nnobio ! Loop over # other surface types
853	!!$ !
854	!!$ IF ( jv .EQ. iice ) THEN
855	!!$ !
856	!!$ ! Land ice is of course a particular case
857	!!$ !
858	!!$ DO ji=1,kjpindex
859	!!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv)
860	!!$ ENDDO
861	!!$ !
862	!!$ ELSE
863	!
864	DO ji=1,kjpindex ! Loop over # pixels - domain size
865
866	vbeta1_add = frac_nobio(ji,jv) * frac_snow_nobio(ji, jv)
867
868	! Limitation of sublimation in case of snow amounts smaller than
869	! the atmospheric demand.
870	speed = MAX(min_wind, wind(ji))
871
872	!! Limitation of sublimation by the snow accumulated on the ground
873	!! A first approximation is obtained with the old values of
874	!! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc)
875	subtest = dt_sechiba * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * &
876	& ( qsatt(ji) - qair(ji) )
877
878	IF ( subtest .GT. min_sechiba ) THEN
879	zrapp = snow_nobio(ji,jv) / subtest
880	IF ( zrapp .LT. un ) THEN
881	vbeta1_add = vbeta1_add * zrapp
882	ENDIF
883	ENDIF
884
885	vbeta1(ji) = vbeta1(ji) + vbeta1_add
886
887	ENDDO ! Loop over # pixels - domain size
888
889	!!$ ENDIF
890
891	ENDDO ! Loop over # other surface types
892
893	IF (printlev>=3) WRITE (numout,*) ' diffuco_snow done '
894
895	END SUBROUTINE diffuco_snow
896
897
898	!! ================================================================================================================================
899	!! SUBROUTINE : diffuco_flood
900	!!
901	!>\BRIEF This routine computes partial beta coefficient : floodplains
902	!!
903	!! DESCRIPTION :
904	!!
905	!! RECENT CHANGE(S): None
906	!!
907	!! MAIN OUTPUT VARIABLE(S) : vbeta5
908	!!
909	!! REFERENCE(S) : None
910	!!
911	!! FLOWCHART : None
912	!! \n
913	!_ ================================================================================================================================
914
915	SUBROUTINE diffuco_flood (kjpindex, qair, qsatt, rau, u, v, q_cdrag, evapot, evapot_corr, &
916	& flood_frac, flood_res, vbeta5)
917
918	! interface description
919	! input scalar
920	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
921	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
922	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity
923	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
924	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
925	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
926	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
927	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_res !! water mass in flood reservoir
928	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! fraction of floodplains
929	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Potential evaporation
930	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr!! Potential evaporation2
931	! output fields
932	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta5 !! Beta for floodplains
933
934	! local declaration
935	REAL(r_std) :: subtest, zrapp, speed
936	INTEGER(i_std) :: ji, jv
937
938	!_ ================================================================================================================================
939	!
940	! beta coefficient for sublimation for floodplains
941	!
942	DO ji=1,kjpindex
943	!
944	IF (evapot(ji) .GT. min_sechiba) THEN
945	vbeta5(ji) = flood_frac(ji) *evapot_corr(ji)/evapot(ji)
946	ELSE
947	vbeta5(ji) = flood_frac(ji)
948	ENDIF
949	!
950	! -- Limitation of evaporation in case of water amounts smaller than
951	! the atmospheric demand.
952
953	!
954	speed = MAX(min_wind, SQRT (u(ji)u(ji) + v(ji)v(ji)))
955	!
956	subtest = dt_sechiba * vbeta5(ji) * speed * q_cdrag(ji) * rau(ji) * &
957	& ( qsatt(ji) - qair(ji) )
958	!
959	IF ( subtest .GT. min_sechiba ) THEN
960	zrapp = flood_res(ji) / subtest
961	IF ( zrapp .LT. un ) THEN
962	vbeta5(ji) = vbeta5(ji) * zrapp
963	ENDIF
964	ENDIF
965	!
966	END DO
967
968	IF (printlev>=3) WRITE (numout,*) ' diffuco_flood done '
969
970	END SUBROUTINE diffuco_flood
971
972
973	!! ================================================================================================================================
974	!! SUBROUTINE : diffuco_inter
975	!!
976	!>\BRIEF This routine computes the partial beta coefficient
977	!! for the interception for each type of vegetation
978	!!
979	!! DESCRIPTION : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax).
980	!! The former is submitted to transpiration only (vbeta3 coefficient, calculated in
981	!! diffuco_trans_co2), while the latter is first submitted to interception loss
982	!! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated
983	!! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test,
984	!! as for snow sublimation. \n
985	!!
986	!! \latexonly
987	!! \input{diffucointer1.tex}
988	!! \endlatexonly
989	!! Calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
990	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
991	!!
992	!! \latexonly
993	!! \input{diffucointer2.tex}
994	!! \endlatexonly
995	!!
996	!! Calculation of $\beta_3$, the canopy transpiration resistance
997	!! \latexonly
998	!! \input{diffucointer3.tex}
999	!! \endlatexonly
1000	!!
1001	!! We here determine the limitation of interception loss by the water stored on the leaf.
1002	!! A first approximation of interception loss is obtained using the old values of
1003	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1004	!! \latexonly
1005	!! \input{diffucointer4.tex}
1006	!! \endlatexonly
1007	!!
1008	!! \latexonly
1009	!! \input{diffucointer5.tex}
1010	!! \endlatexonly
1011	!!
1012	!! \latexonly
1013	!! \input{diffucointer6.tex}
1014	!! \endlatexonly
1015	!!
1016	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1017	!! 'zqsvegrap'.
1018	!! \latexonly
1019	!! \input{diffucointer7.tex}
1020	!! \endlatexonly
1021	!!
1022	!! RECENT CHANGE(S): None
1023	!!
1024	!! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23
1025	!!
1026	!! REFERENCE(S) :
1027	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1028	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1029	!! Circulation Model. Journal of Climate, 6, pp. 248-273
1030	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1031	!! sur le cycle de l'eau, PhD Thesis, available from:
1032	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1033	!! - Perrier, A, 1975. Etude physique de l'Ã©vaporation dans les conditions naturelles. Annales
1034	!! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243
1035	!!
1036	!! FLOWCHART : None
1037	!! \n
1038	!_ ================================================================================================================================
1039
1040	SUBROUTINE diffuco_inter (kjpindex, qair, qsatt, rau, u, v, q_cdrag, humrel, veget, &
1041	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
1042
1043	!! 0 Variable and parameter declaration
1044
1045	!! 0.1 Input variables
1046
1047	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1048	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
1049	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
1050	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
1051	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
1052	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
1053	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1054	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
1055	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type (fraction)
1056	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
1057	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
1058	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! architectural resistance (s m^{-1})
1059
1060	!! 0.2 Output variables
1061
1062	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
1063	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will
1064	!! transpire (-)
1065
1066	!! 0.4 Local variables
1067
1068	INTEGER(i_std) :: ji, jv !! (-), (-)
1069	REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed !!
1070	!_ ================================================================================================================================
1071
1072	!! 1. Initialize
1073
1074	vbeta2(:,:) = zero
1075	vbeta23(:,:) = zero
1076
1077	!! 2. The beta coefficient for interception by vegetation.
1078
1079	DO jv = 2,nvm
1080
1081	DO ji=1,kjpindex
1082
1083	IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN
1084
1085	zqsvegrap = zero
1086	IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN
1087
1088	!! \latexonly
1089	!! \input{diffucointer1.tex}
1090	!! \endlatexonly
1091	!!
1092	!! We calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1093	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1094	!!
1095	zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
1096	END IF
1097
1098	!! \latexonly
1099	!! \input{diffucointer2.tex}
1100	!! \endlatexonly
1101	speed = MAX(min_wind, wind(ji))
1102
1103	!! Calculation of $\beta_3$, the canopy transpiration resistance
1104	!! \latexonly
1105	!! \input{diffucointer3.tex}
1106	!! \endlatexonly
1107	vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv)))
1108
1109	!! We here determine the limitation of interception loss by the water stored on the leaf.
1110	!! A first approximation of interception loss is obtained using the old values of
1111	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1112	!! \latexonly
1113	!! \input{diffucointer4.tex}
1114	!! \endlatexonly
1115	ziltest = dt_sechiba * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * &
1116	& ( qsatt(ji) - qair(ji) )
1117
1118	IF ( ziltest .GT. min_sechiba ) THEN
1119
1120	!! \latexonly
1121	!! \input{diffucointer5.tex}
1122	!! \endlatexonly
1123	zrapp = qsintveg(ji,jv) / ziltest
1124	IF ( zrapp .LT. un ) THEN
1125
1126	!! \latexonly
1127	!! \input{diffucointer6.tex}
1128	!! \endlatexonly
1129	!!
1130	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1131	!! 'zqsvegrap'.
1132	IF ( humrel(ji,jv) >= min_sechiba ) THEN
1133	vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, zero)
1134	ELSE
1135	! We don't want transpiration when the soil cannot deliver it
1136	vbeta23(ji,jv) = zero
1137	ENDIF
1138
1139	!! \latexonly
1140	!! \input{diffucointer7.tex}
1141	!! \endlatexonly
1142	vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp
1143	ENDIF
1144	ENDIF
1145	END IF
1146	! ! Autre formulation possible pour l'evaporation permettant une transpiration sur tout le feuillage
1147	! !commenter si formulation Nathalie sinon Tristan
1148	! speed = MAX(min_wind, wind(ji))
1149	!
1150	! vbeta23(ji,jv) = MAX(zero, veget(ji,jv) * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv))) - vbeta2(ji,jv))
1151
1152	END DO
1153
1154	END DO
1155
1156	IF (printlev>=3) WRITE (numout,*) ' diffuco_inter done '
1157
1158	END SUBROUTINE diffuco_inter
1159
1160
1161	!! ==============================================================================================================================
1162	!! SUBROUTINE : diffuco_bare
1163	!!
1164	!>\BRIEF This routine computes the partial beta coefficient corresponding to
1165	!! bare soil
1166	!!
1167	!! DESCRIPTION : Bare soil evaporation is submitted to a maximum possible flow (evap_bare_lim)
1168	!!
1169	!! Calculation of wind speed
1170	!! \latexonly
1171	!! \input{diffucobare1.tex}
1172	!! \endlatexonly
1173	!!
1174	!! The calculation of $\beta_4$
1175	!! \latexonly
1176	!! \input{diffucobare2.tex}
1177	!! \endlatexonly
1178	!!
1179	!! RECENT CHANGE(S): None
1180	!!
1181	!! MAIN OUTPUT VARIABLE(S): ::vbeta4
1182	!!
1183	!! REFERENCE(S) :
1184	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1185	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1186	!! Circulation Model. Journal of Climate, 6, pp.248-273
1187	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1188	!! sur le cycle de l'eau, PhD Thesis, available from:
1189	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1190	!!
1191	!! FLOWCHART : None
1192	!! \n
1193	!_ ================================================================================================================================
1194
1195	SUBROUTINE diffuco_bare (kjpindex, tot_bare_soil, veget_max, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)
1196
1197	!! 0. Variable and parameter declaration
1198
1199	!! 0.1 Input variables
1200	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1201	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
1202	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
1203	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta2 !! Beta for Interception
1204	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta3 !! Beta for Transpiration
1205
1206	!! 0.2 Output variables
1207	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
1208
1209	!! 0.3 Modified variables
1210	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! limiting factor for bare soil evaporation
1211	!! when the 11-layer hydrology is used (-)
1212	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! limiting factor for bare soil evaporation
1213	!! when the 11-layer hydrology is used (-)
1214	!! 0.4 Local variables
1215	INTEGER(i_std) :: ji
1216	REAL(r_std), DIMENSION(kjpindex) :: vegtot
1217
1218	!_ ================================================================================================================================
1219
1220	!! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil
1221
1222	! To use the new version of hydrol_split_soil and ensure the water conservation, we must have throughout sechiba at all time:
1223	! a) evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)soiltile(ji,:)vegtot(ji))
1224	! b) all the terms (vbeta4, evap_bare_lim, evap_bare_lim_ns) =0 if evap_bare_lim(ji) LE min_sechiba
1225	! This must also be kept true in diffuco_comb
1226
1227	DO ji = 1, kjpindex
1228
1229	! The limitation by 1-beta2-beta3 is due to the fact that evaporation under vegetation is possible
1230	!! \latexonly
1231	!! \input{diffucobare3.tex}
1232	!! \endlatexonly
1233
1234	IF ( (evap_bare_lim(ji) .GT. min_sechiba) .AND. &
1235	! in this case we can't have vegtot LE min_sechina, cf hydrol_soil
1236	(un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) .GT. min_sechiba) ) THEN
1237	! eventually, none of the left-hand term is close to zero
1238
1239	vegtot(ji) = SUM(veget_max(ji,:))
1240	IF (evap_bare_lim(ji) < (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)))) THEN
1241	! Standard case
1242	vbeta4(ji) = evap_bare_lim(ji)
1243	ELSE
1244	vbeta4(ji) = un - SUM(vbeta2(ji,:)+vbeta3(ji,:))
1245
1246	! We now have to redefine evap_bare_lim & evap_bare_lim_ns
1247	IF (evap_bare_lim(ji) .GT. min_sechiba) THEN
1248	evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji)
1249	ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles
1250	evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji)
1251	ENDIF
1252
1253	evap_bare_lim(ji) = vbeta4(ji)
1254	! consistent with evap_bare_lim(ji) =
1255	! SUM(evap_bare_lim_ns(ji,:)soiltile(ji,:)vegtot(ji))
1256	! as SUM(soiltile(ji,:)) = 1
1257	END IF
1258
1259	ELSE ! instead of having very small vbeta4, we set everything to zero
1260	vbeta4(ji) = zero
1261	evap_bare_lim(ji) = zero
1262	evap_bare_lim_ns(ji,:) = zero
1263	ENDIF
1264
1265	END DO
1266
1267	IF (printlev>=3) WRITE (numout,*) ' diffuco_bare done '
1268
1269	END SUBROUTINE diffuco_bare
1270
1271
1272	!! ==============================================================================================================================
1273	!! SUBROUTINE : diffuco_trans_co2
1274	!!
1275	!>\BRIEF This subroutine computes carbon assimilation and stomatal
1276	!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987).
1277	!!
1278	!! DESCRIPTION :\n
1279	!! *** General:\n
1280	!! The equations are different depending on the photosynthesis mode (C3 versus C4).\n
1281	!! Assimilation and conductance are computed over 20 levels of LAI and then
1282	!! integrated at the canopy level.\n
1283	!! This routine also computes partial beta coefficient: transpiration for each
1284	!! type of vegetation.\n
1285	!! There is a main loop on the PFTs, then inner loops on the points where
1286	!! assimilation has to be calculated.\n
1287	!! This subroutine is called at each sechiba time step by sechiba_main.\n
1288	!! *** Details:
1289	!! - Integration at the canopy level\n
1290	!! \latexonly
1291	!! \input{diffuco_trans_co2_1.1.tex}
1292	!! \endlatexonly
1293	!! - Light''s extinction \n
1294	!! The available light follows a simple Beer extinction law.
1295	!! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.\n
1296	!! \latexonly
1297	!! \input{diffuco_trans_co2_1.2.tex}
1298	!! \endlatexonly
1299	!! - Estimation of relative humidity of air (for calculation of the stomatal conductance)\n
1300	!! \latexonly
1301	!! \input{diffuco_trans_co2_1.3.tex}
1302	!! \endlatexonly
1303	!! - Calculation of the water limitation factor\n
1304	!! \latexonly
1305	!! \input{diffuco_trans_co2_2.1.tex}
1306	!! \endlatexonly
1307	!! - Calculation of temperature dependent parameters for C4 plants\n
1308	!! \latexonly
1309	!! \input{diffuco_trans_co2_2.2.tex}
1310	!! \endlatexonly
1311	!! - Calculation of temperature dependent parameters for C3 plants\n
1312	!! \latexonly
1313	!! \input{diffuco_trans_co2_2.3.tex}
1314	!! \endlatexonly
1315	!! - Vmax scaling\n
1316	!! Vmax is scaled into the canopy due to reduction of nitrogen
1317	!! (Johnson and Thornley,1984).\n
1318	!! \latexonly
1319	!! \input{diffuco_trans_co2_2.4.1.tex}
1320	!! \endlatexonly
1321	!! - Assimilation for C4 plants (Collatz et al., 1992)\n
1322	!! \latexonly
1323	!! \input{diffuco_trans_co2_2.4.2.tex}
1324	!! \endlatexonly
1325	!! - Assimilation for C3 plants (Farqhuar et al., 1980)\n
1326	!! \latexonly
1327	!! \input{diffuco_trans_co2_2.4.3.tex}
1328	!! \endlatexonly
1329	!! - Estimation of the stomatal conductance (Ball et al., 1987)\n
1330	!! \latexonly
1331	!! \input{diffuco_trans_co2_2.4.4.tex}
1332	!! \endlatexonly
1333	!!
1334	!! RECENT CHANGE(S):
1335	!!
1336	!! 2018/11 replaced t2m and q2m by temp_air and qair; prognostic variables on first atmospheric model layer
1337	!!
1338	!! 2006/06 N. de Noblet
1339	!! - addition of q2m and t2m as input parameters for the
1340	!! calculation of Rveget
1341	!! - introduction of vbeta23
1342	!!
1343	!! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular
1344	!! concentration
1345	!!
1346	!! REFERENCE(S) :
1347	!! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal
1348	!! conductance and its contribution to the control of photosynthesis under
1349	!! different environmental conditions, Prog. Photosynthesis, 4, 221â 224.
1350	!! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis
1351	!! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol.,
1352	!! 19, 519â538.
1353	!! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of
1354	!! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78â90.
1355	!! - Johnson, I. R., and J. Thornley (1984), A model of instantaneous and daily
1356	!! canopy photosynthesis, J Theor. Biol., 107, 531 545
1357	!! - McMurtrie, R.E., Rook, D.A. and Kelliher, F.M., 1990. Modelling the yield of Pinus radiata on a
1358	!! site limited by water and nitrogen. For. Ecol. Manage., 30: 381-413
1359	!! - Bounoua, L., Hall, F. G., Sellers, P. J., Kumar, A., Collatz, G. J., Tucker, C. J., and Imhoff, M. L. (2010), Quantifying the
1360	!! negative feedback of vegetation to greenhouse warming: A modeling approach, Geophysical Research Letters, 37, Artn L23701,
1361	!! Doi 10.1029/2010gl045338
1362	!! - Bounoua, L., Collatz, G. J., Sellers, P. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I.,
1363	!! Tucker, C. J., Field, C. B., and Jensen, T. G. (1999), Interactions between vegetation and climate: Radiative and physiological
1364	!! effects of doubled atmospheric co2, Journal of Climate, 12, 309-324, Doi 10.1175/1520-0442(1999)012<0309:Ibvacr>2.0.Co;2
1365	!! - Sellers, P. J., Bounoua, L., Collatz, G. J., Randall, D. A., Dazlich, D. A., Los, S. O., Berry, J. A., Fung, I.,
1366	!! Tucker, C. J., Field, C. B., and Jensen, T. G. (1996), Comparison of radiative and physiological effects of doubled atmospheric
1367	!! co2 on climate, Science, 271, 1402-1406, DOI 10.1126/science.271.5254.1402
1368	!! - Lewis, J. D., Ward, J. K., and Tissue, D. T. (2010), Phosphorus supply drives nonlinear responses of cottonwood
1369	!! (populus deltoides) to increases in co2 concentration from glacial to future concentrations, New Phytologist, 187, 438-448,
1370	!! DOI 10.1111/j.1469-8137.2010.03307.x
1371	!! - Kattge, J., Knorr, W., Raddatz, T., and Wirth, C. (2009), Quantifying photosynthetic capacity and its relationship to leaf
1372	!! nitrogen content for global-scale terrestrial biosphere models, Global Change Biology, 15, 976-991,
1373	!! DOI 10.1111/j.1365-2486.2008.01744.x
1374	!!
1375	!! FLOWCHART : None
1376	!! \n
1377	!_ ================================================================================================================================
1378
1379	SUBROUTINE diffuco_trans_co2 (kjpindex, swdown, pb, qsurf, qair, temp_air, temp_growth, rau, u, v, q_cdrag, humrel, &
1380	assim_param, Ca, &
1381	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, vbeta3pot, rveget, rstruct, &
1382	cimean, gsmean, gpp, &
1383	vbeta23, hist_id, indexveg, indexlai, index, kjit, cim)
1384
1385	!
1386	!! 0. Variable and parameter declaration
1387	!
1388
1389	!
1390	!! 0.1 Input variables
1391	!
1392	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1393	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Downwelling short wave flux
1394	!! @tex ($W m^{-2}$) @endtex
1395	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
1396	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface specific humidity
1397	!! @tex ($kg kg^{-1}$) @endtex
1398	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Specific humidity at first atmospheric model layer
1399	!! @tex ($kg kg^{-1}$) @endtex
1400	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature at first atmospheric model layer (K)
1401	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_growth !! Growth temperature (Â°C) - Is equal to t2m_month
1402	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! air density @tex ($kg m^{-3}$) @endtex
1403	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
1404	!! @tex ($m s^{-1}$) @endtex
1405	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
1406	!! @tex ($m s^{-1}$) @endtex
1407	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
1408	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! min+max+opt temps (K), vcmax, vjmax for
1409	!! photosynthesis
1410	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1411	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: Ca !! CO2 concentration inside the canopy
1412	!! @tex ($\mu mol mol^{-1}$) @endtex
1413	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (0-1,unitless)
1414	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Coverage fraction of vegetation for each PFT
1415	!! depending on LAI (0-1, unitless)
1416	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum vegetation fraction of each PFT inside
1417	!! the grid box (0-1, unitless)
1418	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index @tex ($m^2 m^{-2}$) @endtex
1419	!! @tex ($m s^{-1}$) @endtex
1420	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
1421	!! @tex ($kg m^{-2}$) @endte
1422	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
1423	!! @tex ($kg m^{-2}$) @endtex
1424	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will
1425	!! transpire (unitless)
1426	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-)
1427	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-)
1428	INTEGER(i_std),DIMENSION (kjpindex*(nlai+1)), INTENT (in) :: indexlai !! Indeces of the points on the 3D map
1429	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
1430	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
1431	!
1432	!! 0.2 Output variables
1433	!
1434	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration (unitless)
1435	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3pot !! Beta for Potential Transpiration
1436	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! stomatal resistance of vegetation
1437	!! @tex ($s m^{-1}$) @endtex
1438	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex
1439	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! mean intercellular CO2 concentration
1440	!! @tex ($\mu mol mol^{-1}$) @endtex
1441	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: gsmean !! mean stomatal conductance to CO2 (umol m-2 s-1)
1442	REAL(r_Std),DIMENSION (kjpindex,nvm), INTENT (out) :: gpp !! Assimilation ((gC m^{-2} dt_sechiba^{-1}), total area)
1443	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cim !! Intercellular CO2 over nlai
1444	!
1445	!! 0.3 Modified variables
1446	!
1447
1448	!
1449	!! 0.4 Local variables
1450	!
1451	REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax !! maximum rate of carboxylation
1452	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1453	INTEGER(i_std) :: ji, jv, jl, limit_photo !! indices (unitless)
1454	REAL(r_std), DIMENSION(kjpindex,nlai+1) :: info_limitphoto
1455	REAL(r_std), DIMENSION(kjpindex,nvm,nlai) :: leaf_ci !! intercellular CO2 concentration (ppm)
1456	REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest !! intercellular CO2 concentration at the lowest
1457	!! LAI level
1458	!! @tex ($\mu mol mol^{-1}$) @endtex
1459	INTEGER(i_std), DIMENSION(kjpindex) :: ilai !! counter for loops on LAI levels (unitless)
1460	REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap !! relative water quantity in the water
1461	!! interception reservoir (0-1,unitless)
1462	REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex
1463	! Assimilation
1464	LOGICAL, DIMENSION(kjpindex) :: assimilate !! where assimilation is to be calculated
1465	!! (unitless)
1466	LOGICAL, DIMENSION(kjpindex) :: calculate !! where assimilation is to be calculated for
1467	!! in the PFTs loop (unitless)
1468	INTEGER(i_std) :: nic,inic,icinic !! counter/indices (unitless)
1469	INTEGER(i_std), DIMENSION(kjpindex) :: index_calc !! index (unitless)
1470	INTEGER(i_std) :: nia,inia,nina,inina,iainia !! counter/indices (unitless)
1471	INTEGER(i_std), DIMENSION(kjpindex) :: index_assi,index_non_assi !! indices (unitless)
1472	REAL(r_std), DIMENSION(kjpindex, nlai+1) :: vc2 !! rate of carboxylation (at a specific LAI level)
1473	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1474	REAL(r_std), DIMENSION(kjpindex, nlai+1) :: vj2 !! rate of Rubisco regeneration (at a specific LAI
1475	!! level) @tex ($\mu mol e- m^{-2} s^{-1}$) @endtex
1476	REAL(r_std), DIMENSION(kjpindex, nlai+1) :: assimi !! assimilation (at a specific LAI level)
1477	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1478	!! (temporary variables)
1479	REAL(r_std), DIMENSION(kjpindex) :: gstop !! stomatal conductance to H2O at topmost level
1480	!! @tex ($m s^{-1}$) @endtex
1481	REAL(r_std), DIMENSION(kjpindex,nlai+1) :: gs !! stomatal conductance to CO2
1482	!! @tex ($\mol m^{-2} s^{-1}$) @endtex
1483	REAL(r_std), DIMENSION(kjpindex,nlai) :: templeafci
1484
1485
1486	REAL(r_std), DIMENSION(kjpindex) :: gamma_star !! CO2 compensation point (ppm)
1487	!! @tex ($\mu mol mol^{-1}$) @endtex
1488
1489	REAL(r_std), DIMENSION(kjpindex) :: air_relhum !! air relative humidity at 2m
1490	!! @tex ($kg kg^{-1}$) @endtex
1491	REAL(r_std), DIMENSION(kjpindex) :: VPD !! Vapor Pressure Deficit (kPa)
1492	REAL(r_std), DIMENSION(kjpindex) :: water_lim !! water limitation factor (0-1,unitless)
1493
1494	REAL(r_std), DIMENSION(kjpindex) :: gstot !! total stomatal conductance to H2O
1495	!! Final unit is
1496	!! @tex ($m s^{-1}$) @endtex
1497	REAL(r_std), DIMENSION(kjpindex) :: assimtot !! total assimilation
1498	!! @tex ($\mu mol CO2 m^{-2} s^{-1}$) @endtex
1499	REAL(r_std), DIMENSION(kjpindex) :: Rdtot !! Total Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1)
1500	REAL(r_std), DIMENSION(kjpindex) :: leaf_gs_top !! leaf stomatal conductance to H2O at topmost level
1501	!! @tex ($\mol H2O m^{-2} s^{-1}$) @endtex
1502	REAL(r_std), DIMENSION(nlai+1) :: laitab !! tabulated LAI steps @tex ($m^2 m^{-2}$) @endtex
1503	REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??)
1504	!! @tex ($g g^{-1}$) @endtex
1505	REAL(r_std), DIMENSION(nvm,nlai) :: light !! fraction of light that gets through upper LAI
1506	!! levels (0-1,unitless)
1507	REAL(r_std), DIMENSION(kjpindex) :: T_Vcmax !! Temperature dependance of Vcmax (unitless)
1508	REAL(r_std), DIMENSION(kjpindex) :: S_Vcmax_acclim_temp !! Entropy term for Vcmax
1509	!! accounting for acclimation to temperature (J K-1 mol-1)
1510	REAL(r_std), DIMENSION(kjpindex) :: T_Jmax !! Temperature dependance of Jmax
1511	REAL(r_std), DIMENSION(kjpindex) :: S_Jmax_acclim_temp !! Entropy term for Jmax
1512	!! accounting for acclimation toxs temperature (J K-1 mol-1)
1513	REAL(r_std), DIMENSION(kjpindex) :: T_gm !! Temperature dependance of gmw
1514	REAL(r_std), DIMENSION(kjpindex) :: T_Rd !! Temperature dependance of Rd (unitless)
1515	REAL(r_std), DIMENSION(kjpindex) :: T_Kmc !! Temperature dependance of KmC (unitless)
1516	REAL(r_std), DIMENSION(kjpindex) :: T_KmO !! Temperature dependance of KmO (unitless)
1517	REAL(r_std), DIMENSION(kjpindex) :: T_Sco !! Temperature dependance of Sco
1518	REAL(r_std), DIMENSION(kjpindex) :: T_gamma_star !! Temperature dependance of gamma_star (unitless)
1519	REAL(r_std), DIMENSION(kjpindex) :: vc !! Maximum rate of Rubisco activity-limited carboxylation (mumol CO2 mâ2 sâ1)
1520	REAL(r_std), DIMENSION(kjpindex) :: vj !! Maximum rate of e- transport under saturated light (mumol CO2 mâ2 sâ1)
1521	REAL(r_std), DIMENSION(kjpindex) :: gm !! Mesophyll diffusion conductance (molCO2 mâ2 sâ1 barâ1)
1522	REAL(r_std), DIMENSION(kjpindex) :: g0var
1523	REAL(r_std), DIMENSION(kjpindex,nlai+1) :: Rd !! Day respiration (respiratory CO2 release other than by photorespiration) (mumol CO2 mâ2 sâ1)
1524	REAL(r_std), DIMENSION(kjpindex) :: Kmc !! MichaelisâMenten constant of Rubisco for CO2 (mubar)
1525	REAL(r_std), DIMENSION(kjpindex) :: KmO !! MichaelisâMenten constant of Rubisco for O2 (mubar)
1526	REAL(r_std), DIMENSION(kjpindex) :: Sco !! Relative CO2 /O2 specificity factor for Rubisco (bar bar-1)
1527	REAL(r_std), DIMENSION(kjpindex) :: gb_co2 !! Boundary-layer conductance (molCO2 mâ2 sâ1 barâ1)
1528	REAL(r_std), DIMENSION(kjpindex) :: gb_h2o !! Boundary-layer conductance (molH2O mâ2 sâ1 barâ1)
1529	REAL(r_std), DIMENSION(kjpindex) :: fvpd !! Factor for describing the effect of leaf-to-air vapour difference on gs (-)
1530	REAL(r_std), DIMENSION(kjpindex) :: low_gamma_star !! Half of the reciprocal of Sc/o (bar bar-1)
1531	REAL(r_std) :: N_Vcmax !! Nitrogen level dependance of Vcmacx and Jmax
1532	REAL(r_std) :: fcyc !! Fraction of electrons at PSI that follow cyclic transport around PSI (-)
1533	REAL(r_std) :: z !! A lumped parameter (see Yin et al. 2009) ( mol mol-1)
1534	REAL(r_std) :: Rm !! Day respiration in the mesophyll (umol CO2 mâ2 sâ1)
1535	REAL(r_std) :: Cs_star !! Cs -based CO2 compensation point in the absence of Rd (ubar)
1536	REAL(r_std), DIMENSION(kjpindex) :: Iabs !! Photon flux density absorbed by leaf photosynthetic pigments (umol photon mâ2 sâ1)
1537	REAL(r_std), DIMENSION(kjpindex) :: Jmax !! Maximum value of J under saturated light (umol eâ mâ2 sâ1)
1538	REAL(r_std), DIMENSION(kjpindex, nlai+1) :: JJ !! Rate of eâ transport (umol eâ mâ2 sâ1)
1539	REAL(r_std) :: J2 !! Rate of all eâ transport through PSII (umol eâ mâ2 sâ1)
1540	REAL(r_std) :: VpJ2 !! eâ transport-limited PEP carboxylation rate (umol CO2 mâ2 sâ1)
1541	REAL(r_std) :: A_1, A_3 !! Lowest First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
1542	REAL(r_std) :: A_1_tmp, A_3_tmp !! Temporary First and third roots of the analytical solution for a general cubic equation (see Appendix A of Yin et al. 2009) (umol CO2 mâ2 sâ1)
1543	REAL(r_std) :: Obs !! Bundle-sheath oxygen partial pressure (ubar)
1544	REAL(r_std), DIMENSION(kjpindex, nlai+1) :: Cc !! Chloroplast CO2 partial pressure (ubar)
1545	REAL(r_std) :: ci_star !! Ci -based CO2 compensation point in the absence of Rd (ubar)
1546	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))
1547	REAL(r_std) :: QQ,UU,PSI,x1,x2,x3 !! Variables used for solving the cubic equation (see Yin et al. (2009))
1548
1549	REAL(r_std) :: cresist !! coefficient for resistances (??)
1550	REAL(r_std), DIMENSION(kjpindex) :: laisum !! when calculating cim over nlai
1551
1552	! @defgroup Photosynthesis Photosynthesis
1553	! @{
1554	! 1. Preliminary calculations\n
1555	!_ ================================================================================================================================
1556
1557	cim(:,:)=zero
1558	leaf_ci(:,:,:) = zero
1559
1560	!
1561	! 1.1 Calculate LAI steps\n
1562	! The integration at the canopy level is done over nlai fixed levels.
1563	!! \latexonly
1564	!! \input{diffuco_trans_co2_1.1.tex}
1565	!! \endlatexonly
1566	! @}
1567	! @codeinc
1568	DO jl = 1, nlai+1
1569	laitab(jl) = laimax(EXP(lai_level_depthREAL(jl-1,r_std))-1.)/(EXP(lai_level_depth*REAL(nlai,r_std))-un)
1570	ENDDO
1571	! @endcodeinc
1572
1573	! @addtogroup Photosynthesis
1574	! @{
1575	!
1576	! 1.2 Calculate light fraction for each LAI step\n
1577	! The available light follows a simple Beer extinction law.
1578	! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.
1579	!! \latexonly
1580	!! \input{diffuco_trans_co2_1.2.tex}
1581	!! \endlatexonly
1582	! @}
1583	! @codeinc
1584	DO jl = 1, nlai
1585	DO jv = 1, nvm
1586	light(jv,jl) = exp( -ext_coeff(jv)*laitab(jl) )
1587	ENDDO
1588	ENDDO
1589	! @endcodeinc
1590	!
1591	! Photosynthesis parameters
1592	!
1593
1594	! Choice of downregulation. Note that downregulation_co2_new excludes
1595	IF (downregulation_co2_new) THEN
1596	! Option used for CMIP6 version from 6.1.11
1597	! A minimum value is used for the CO2 concentration(Ca)
1598	! For low CO2 concentration values(under downregulation_co2_minimum) the
1599	! parametrization allows a behavior in the same way as for the preindustral period.
1600	DO jv= 1, nvm
1601	vcmax(:,jv) = assim_param(:,jv,ivcmax)(un-downregulation_co2_coeff_new(jv) &
1602	(MAX(Ca(:),downregulation_co2_minimum)-downregulation_co2_baselevel) / &
1603	(MAX(Ca(:),downregulation_co2_minimum)+20.))
1604	ENDDO
1605	ELSE IF (downregulation_co2) THEN
1606	! Option used for CMIP6 version 6.1.0 up to 6.1.10
1607	DO jv= 1, nvm
1608	vcmax(:,jv) = assim_param(:,jv,ivcmax)(un-downregulation_co2_coeff(jv) &
1609	log(Ca(:)/downregulation_co2_baselevel))
1610	ENDDO
1611	ELSE
1612	vcmax(:,:) = assim_param(:,:,ivcmax)
1613	ENDIF
1614
1615	! DO jv = 1, nvm
1616	! vcmax(:,:) = Vcmax25(jv)
1617	! ENDDO
1618
1619	! @addtogroup Photosynthesis
1620	! @{
1621	!
1622	! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance).\n
1623	!! \latexonly
1624	!! \input{diffuco_trans_co2_1.3.tex}
1625	!! \endlatexonly
1626	! @}
1627	!
1628
1629	CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
1630	air_relhum(:) = &
1631	( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) ) / &
1632	( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) )
1633
1634
1635	VPD(:) = ( qsatt(:)pb(:) / (Tetens_1+qsatt(:)Tetens_2 ) ) &
1636	- ( qair(:) * pb(:) / (Tetens_1+qair(:)* Tetens_2) )
1637	! VPD is needed in kPa
1638	VPD(:) = VPD(:)/10.
1639
1640	!
1641	! 2. beta coefficient for vegetation transpiration
1642	!
1643	rstruct(:,1) = rstruct_const(1)
1644	rveget(:,:) = undef_sechiba
1645	!
1646	vbeta3(:,:) = zero
1647	vbeta3pot(:,:) = zero
1648	gsmean(:,:) = zero
1649	gpp(:,:) = zero
1650	!
1651	cimean(:,1) = Ca(:)
1652	!
1653	! @addtogroup Photosynthesis
1654	! @{
1655	! 2. Loop over vegetation types\n
1656	! @}
1657	!
1658	DO jv = 2,nvm
1659	gamma_star(:)=zero
1660	Kmo(:)=zero
1661	Kmc(:)=zero
1662	gm(:)=zero
1663	g0var(:) =zero
1664
1665	Cc(:,:)=zero
1666	Vc2(:,:)=zero
1667	JJ(:,:)=zero
1668	info_limitphoto(:,:)=zero
1669	gs(:,:)=zero
1670	templeafci(:,:)=zero
1671	assimi(:,:)=zero
1672	Rd(:,:)=zero
1673
1674	!
1675	! @addtogroup Photosynthesis
1676	! @{
1677	!
1678	! 2.1 Initializations\n
1679	!! \latexonly
1680	!! \input{diffuco_trans_co2_2.1.tex}
1681	!! \endlatexonly
1682	! @}
1683	!
1684	! beta coefficient for vegetation transpiration
1685	!
1686	rstruct(:,jv) = rstruct_const(jv)
1687	cimean(:,jv) = Ca(:)
1688	!
1689	!! mask that contains points where there is photosynthesis
1690	!! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points.
1691	!! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not
1692	!! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation,
1693	!! temperature and relative humidity).
1694	!! For the points where assimilation is not calculated, variables are initialized to specific values.
1695	!! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
1696	!! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes
1697	!! the nina points with no assimilation.
1698	nia=0
1699	nina=0
1700	!
1701	DO ji=1,kjpindex
1702	!
1703	IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. &
1704	( veget_max(ji,jv) .GT. min_sechiba ) ) THEN
1705
1706	IF ( ( veget(ji,jv) .GT. min_sechiba ) .AND. &
1707	( swdown(ji) .GT. min_sechiba ) .AND. &
1708	( humrel(ji,jv) .GT. min_sechiba) .AND. &
1709	( temp_growth(ji) .GT. tphoto_min(jv) ) .AND. &
1710	( temp_growth(ji) .LT. tphoto_max(jv) ) ) THEN
1711	!
1712	assimilate(ji) = .TRUE.
1713	nia=nia+1
1714	index_assi(nia)=ji
1715	!
1716	ELSE
1717	!
1718	assimilate(ji) = .FALSE.
1719	nina=nina+1
1720	index_non_assi(nina)=ji
1721	!
1722	ENDIF
1723	ELSE
1724	!
1725	assimilate(ji) = .FALSE.
1726	nina=nina+1
1727	index_non_assi(nina)=ji
1728	!
1729	ENDIF
1730	!
1731
1732	ENDDO
1733	!
1734
1735	gstot(:) = zero
1736	gstop(:) = zero
1737	assimtot(:) = zero
1738	Rdtot(:)=zero
1739	leaf_gs_top(:) = zero
1740	!
1741	zqsvegrap(:) = zero
1742	WHERE (qsintmax(:,jv) .GT. min_sechiba)
1743	!! relative water quantity in the water interception reservoir
1744	zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
1745	ENDWHERE
1746	!
1747	!! Calculates the water limitation factor.
1748	water_lim(:) = humrel(:,jv)
1749
1750	! give a default value of ci for all pixel that do not assimilate
1751	DO jl=1,nlai
1752	DO inina=1,nina
1753	leaf_ci(index_non_assi(inina),jv,jl) = Ca(index_non_assi(inina))
1754	ENDDO
1755	ENDDO
1756	!
1757	ilai(:) = 1
1758	!
1759	! Here is the calculation of assimilation and stomatal conductance
1760	! based on the work of Farquahr, von Caemmerer and Berry (FvCB model)
1761	! as described in Yin et al. 2009
1762	! Yin et al. developed a extended version of the FvCB model for C4 plants
1763	! and proposed an analytical solution for both photosynthesis pathways (C3 and C4)
1764	! Photosynthetic parameters used are those reported in Yin et al.
1765	! Except For Vcmax25, relationships between Vcmax25 and Jmax25 for which we use
1766	! Medlyn et al. (2002) and Kattge & Knorr (2007)
1767	! Because these 2 references do not consider mesophyll conductance, we neglect this term
1768	! in the formulations developed by Yin et al.
1769	! Consequently, gm (the mesophyll conductance) tends to the infinite
1770	! This is of importance because as stated by Kattge & Knorr and Medlyn et al.,
1771	! values of Vcmax and Jmax derived with different model parametrizations are not
1772	! directly comparable and the published values of Vcmax and Jmax had to be standardized
1773	! to one consistent formulation and parametrization
1774
1775	! See eq. 6 of Yin et al. (2009)
1776	! Parametrization of Medlyn et al. (2002) - from Bernacchi et al. (2001)
1777	T_KmC(:) = Arrhenius(kjpindex,temp_air,298.,E_KmC(jv))
1778	T_KmO(:) = Arrhenius(kjpindex,temp_air,298.,E_KmO(jv))
1779	T_Sco(:) = Arrhenius(kjpindex,temp_air,298.,E_Sco(jv))
1780	T_gamma_star(:) = Arrhenius(kjpindex,temp_air,298.,E_gamma_star(jv))
1781
1782
1783	! Parametrization of Yin et al. (2009) - from Bernacchi et al. (2001)
1784	T_Rd(:) = Arrhenius(kjpindex,temp_air,298.,E_Rd(jv))
1785
1786
1787	! For C3 plants, we assume that the Entropy term for Vcmax and Jmax
1788	! acclimates to temperature as shown by Kattge & Knorr (2007) - Eq. 9 and 10
1789	! and that Jmax and Vcmax respond to temperature following a modified Arrhenius function
1790	! (with a decrease of these parameters for high temperature) as in Medlyn et al. (2002)
1791	! and Kattge & Knorr (2007).
1792	! In Yin et al. (2009), temperature dependance to Vcmax is based only on a Arrhenius function
1793	! Concerning this apparent unconsistency, have a look to the section 'Limitation of
1794	! Photosynthesis by gm' of Bernacchi (2002) that may provide an explanation
1795
1796	! Growth temperature tested by Kattge & Knorr range from 11 to 35Â°C
1797	! So, we limit the relationship between these lower and upper limits
1798	S_Jmax_acclim_temp(:) = aSJ(jv) + bSJ(jv) * MAX(11., MIN(temp_growth(:),35.))
1799	T_Jmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Jmax(jv),D_Jmax(jv),S_Jmax_acclim_temp)
1800
1801	S_Vcmax_acclim_temp(:) = aSV(jv) + bSV(jv) * MAX(11., MIN(temp_growth(:),35.))
1802	T_Vcmax(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_Vcmax(jv),D_Vcmax(jv),S_Vcmax_acclim_temp)
1803
1804
1805
1806	vc(:) = vcmax(:,jv) * T_Vcmax(:)
1807	! As shown by Kattge & Knorr (2007), we make use
1808	! of Jmax25/Vcmax25 ratio (rJV) that acclimates to temperature for C3 plants
1809	! rJV is written as a function of the growth temperature
1810	! rJV = arJV + brJV * T_month
1811	! See eq. 10 of Kattge & Knorr (2007)
1812	! and Table 3 for Values of arJV anf brJV
1813	! Growth temperature is monthly temperature (expressed in Â°C) - See first paragraph of
1814	! section Methods/Data of Kattge & Knorr
1815	vj(:) = ( arJV(jv) + brJV(jv) * MAX(11., MIN(temp_growth(:),35.)) ) * vcmax(:,jv) * T_Jmax(:)
1816
1817	T_gm(:) = Arrhenius_modified(kjpindex,temp_air,298.,E_gm(jv),D_gm(jv),S_gm(jv))
1818	gm(:) = gm25(jv) * T_gm(:) * MAX(1-stress_gm(jv), water_lim(:))
1819
1820	g0var(:) = g0(jv)* MAX(1-stress_gs(jv), water_lim(:))
1821	! @endcodeinc
1822	!
1823	KmC(:)=KmC25(jv)*T_KmC(:)
1824	KmO(:)=KmO25(jv)*T_KmO(:)
1825	Sco(:)=Sco25(jv)*T_sco(:)
1826	gamma_star(:) = gamma_star25(jv)*T_gamma_star(:)
1827
1828
1829
1830	! low_gamma_star is defined by Yin et al. (2009)
1831	! as the half of the reciprocal of Sco - See Table 2
1832	low_gamma_star(:) = 0.5 / Sco(:)
1833
1834	! VPD expressed in kPa
1835	! Note : MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) is always between 0-1 not including 0 and 1
1836	fvpd(:) = 1. / ( 1. / MIN(1.-min_sechiba,MAX(min_sechiba,(a1(jv) - b1(jv) * VPD(:)))) - 1. ) &
1837	* MAX(1-stress_gs(jv), water_lim(:))
1838
1839	! leaf boundary layer conductance
1840	! conversion from a conductance in (m s-1) to (mol H2O m-2 s-1)
1841	! from Pearcy et al. (1991, see below)
1842	gb_h2o(:) = gb_ref * 44.6 * (tp_00/temp_air(:)) * (pb(:)/pb_std)
1843
1844	! conversion from (mol H2O m-2 s-1) to (mol CO2 m-2 s-1)
1845	gb_co2(:) = gb_h2o(:) / ratio_H2O_to_CO2
1846
1847	!
1848	! @addtogroup Photosynthesis
1849	! @{
1850	!
1851	! 2.4 Loop over LAI discretized levels to estimate assimilation and conductance\n
1852	! @}
1853	!
1854	!! The calculate(kjpindex) array is of type logical to indicate wether we have to sum over this LAI fixed level (the LAI of
1855	!! the point for the PFT is lower or equal to the LAI level value). The number of such points is incremented in nic and the
1856	!! corresponding point is indexed in the index_calc array.
1857	JJ(:,:)=zero
1858	vc2(:,:)=zero
1859	vj2(:,:)=zero
1860	Cc(:,:)=zero
1861	gs(:,:)=zero
1862	assimi(:,:)=zero
1863	Rd(:,:)=zero
1864
1865	DO jl = 1, nlai
1866	!
1867	nic=0
1868	calculate(:) = .FALSE.
1869	!
1870	IF (nia .GT. 0) then
1871	DO inia=1,nia
1872	calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) )
1873	IF ( calculate(index_assi(inia)) ) THEN
1874	nic=nic+1
1875	index_calc(nic)=index_assi(inia)
1876	ENDIF
1877	ENDDO
1878	ENDIF
1879	!
1880	! @addtogroup Photosynthesis
1881	! @{
1882	!
1883	! 2.4.1 Vmax is scaled into the canopy due to reduction of nitrogen
1884	!! (Johnson and Thornley,1984).\n
1885	!! \latexonly
1886	!! \input{diffuco_trans_co2_2.4.1.tex}
1887	!! \endlatexonly
1888	! @}
1889	!
1890	N_Vcmax = ( un - .7_r_std * ( un - light(jv,jl) ) )
1891	!
1892
1893	vc2(:,jl) = vc(:) * N_Vcmax * MAX(1-stress_vcmax(jv), water_lim(:))
1894	vj2(:,jl) = vj(:) * N_Vcmax * MAX(1-stress_vcmax(jv), water_lim(:))
1895
1896	! see Comment in legend of Fig. 6 of Yin et al. (2009)
1897	! Rd25 is assumed to equal 0.01 Vcmax25
1898	Rd(:,jl) = vcmax(:,jv) * N_Vcmax * 0.01 * T_Rd(:) * MAX(1-stress_vcmax(jv), water_lim(:))
1899
1900	Iabs(:)=swdown(:)W_to_molRG_to_PARext_coeff(jv)light(jv,jl)
1901
1902	! eq. 4 of Yin et al (2009)
1903	Jmax(:)=vj2(:,jl)
1904	JJ(:,jl) = ( alpha_LL(jv) * Iabs(:) + Jmax(:) - sqrt((alpha_LL(jv) * Iabs(:) + Jmax(:) )**2. &
1905	- 4 * theta(jv) * Jmax(:) * alpha_LL(jv) * Iabs(:)) ) &
1906	/ ( 2 * theta(jv))
1907
1908	!
1909	IF ( is_c4(jv) ) THEN
1910	!
1911	! @addtogroup Photosynthesis
1912	! @{
1913	!
1914	! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n
1915	!! \latexonly
1916	!! \input{diffuco_trans_co2_2.4.2.tex}
1917	!! \endlatexonly
1918	! @}
1919	!
1920	!
1921	!
1922	IF (nic .GT. 0) THEN
1923	DO inic=1,nic
1924
1925	! Analytical resolution of the Assimilation based Yin et al. (2009)
1926	icinic=index_calc(inic)
1927
1928	! Eq. 28 of Yin et al. (2009)
1929	fcyc= 1. - ( 4.(1.-fpsir(jv))(1.+fQ(jv)) + 3.h_protons(jv)fpseudo(jv) ) / &
1930	( 3.h_protons(jv) - 4.(1.-fpsir(jv)))
1931
1932	! See paragraph after eq. (20b) of Yin et al.
1933	Rm=Rd(icinic,jl)/2.
1934
1935	! We assume that cs_star equals ci_star (see Comment in legend of Fig. 6 of Yin et al. (2009)
1936	! Equation 26 of Yin et al. (2009)
1937	Cs_star = (gbs(jv) * low_gamma_star(icinic) * Oi - &
1938	( 1. + low_gamma_star(icinic) * alpha(jv) / 0.047) * Rd(icinic,jl) + Rm ) &
1939	/ ( gbs(jv) + kp(jv) )
1940
1941	! eq. 11 of Yin et al (2009)
1942	J2 = JJ(icinic,jl) / ( 1. - fpseudo(jv) / ( 1. - fcyc ) )
1943
1944	! Equation right after eq. (20d) of Yin et al. (2009)
1945	z = ( 2. + fQ(jv) - fcyc ) / ( h_protons(jv) * (1. - fcyc ))
1946
1947	VpJ2 = fpsir(jv) * J2 * z / 2.
1948
1949	A_3=9999.
1950
1951	! See eq. right after eq. 18 of Yin et al. (2009)
1952	DO limit_photo=1,2
1953	! Is Vc limiting the Assimilation
1954	IF ( limit_photo .EQ. 1 ) THEN
1955	a = 1. + kp(jv) / gbs(jv)
1956	b = 0.
1957	x1 = Vc2(icinic,jl)
1958	x2 = KmC(icinic)/KmO(icinic)
1959	x3 = KmC(icinic)
1960	! Is J limiting the Assimilation
1961	ELSE
1962	a = 1.
1963	b = VpJ2
1964	x1 = (1.- fpsir(jv)) * J2 * z / 3.
1965	x2 = 7. * low_gamma_star(icinic) / 3.
1966	x3 = 0.
1967	ENDIF
1968
1969	m=fvpd(icinic)-g0var(icinic)/gb_co2(icinic)
1970	d=g0var(icinic)(Ca(icinic)-Cs_star) + fvpd(icinic)Rd(icinic,jl)
1971	f=(b-Rm-low_gamma_star(icinic)Oigbs(jv))x1d + agbs(jv)x1Ca(icinic)d
1972	j=(b-Rm+gbs(jv)x3 + x2gbs(jv)Oi)m + (alpha(jv)x2/0.047-1.)d &
1973	+ agbs(jv)(Ca(icinic)*m - d/gb_co2(icinic) - (Ca(icinic) - Cs_star ))
1974
1975	g=(b-Rm-low_gamma_star(icinic)Oigbs(jv))x1m - (alpha(jv)low_gamma_star(icinic)/0.047+1.)x1*d &
1976	+ agbs(jv)x1(Ca(icinic)m - d/gb_co2(icinic) - (Ca(icinic)-Cs_star ))
1977
1978	h=-((alpha(jv)low_gamma_star(icinic)/0.047+1.)x1m + (agbs(jv)x1(m-1.))/gb_co2(icinic) )
1979	i= ( b-Rm + gbs(jv)x3 + x2gbs(jv)Oi )d + agbs(jv)Ca(icinic)*d
1980	l= ( alpha(jv)x2/0.047 - 1.)m - (agbs(jv)(m-1.))/gb_co2(icinic)
1981
1982	p = (j-(h-l*Rd(icinic,jl))) / l
1983	q = (i+j*Rd(icinic,jl)-g) / l
1984	r = -(f-i*Rd(icinic,jl)) / l
1985
1986	! See Yin et al. (2009) and Baldocchi (1994)
1987	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
1988	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
1989
1990	IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN
1991	PSI = ACOS(UU/(QQ**1.5_r_std))
1992	A_3_tmp = -2._r_std * SQRT(QQ) * COS(( PSI + 4._r_std * PI)/3._r_std ) - p / 3._r_std
1993	IF (( A_3_tmp .LT. A_3 )) THEN
1994	A_3 = A_3_tmp
1995	info_limitphoto(icinic,jl)=2.
1996	ELSE
1997	! In case, J is not limiting the assimilation
1998	! we have to re-initialise a, b, x1, x2 and x3 values
1999	! in agreement with a Vc-limited assimilation
2000	a = 1. + kp(jv) / gbs(jv)
2001	b = 0.
2002	x1 = Vc2(icinic,jl)
2003	x2 = KmC(icinic)/KmO(icinic)
2004	x3 = KmC(icinic)
2005	info_limitphoto(icinic,jl)=1.
2006	ENDIF
2007	ENDIF
2008
2009	IF ( ( A_3 .EQ. 9999. ) .OR. ( A_3 .LT. (-Rd(icinic,jl)) ) ) THEN
2010	IF ( printlev>=4 ) THEN
2011	WRITE(numout,*) 'We have a problem in diffuco_trans_co2 for A_3'
2012	WRITE(numout,*) 'no real positive solution found for pft:',jv
2013	WRITE(numout,*) 'temp_air:',temp_air(icinic)
2014	WRITE(numout,*) 'vpd:',vpd(icinic)
2015	END IF
2016	A_3 = -Rd(icinic,jl)
2017	ENDIF
2018	assimi(icinic,jl) = A_3
2019
2020	IF ( ABS( assimi(icinic,jl) + Rd(icinic,jl) ) .LT. min_sechiba ) THEN
2021	gs(icinic,jl) = g0var(icinic)
2022	!leaf_ci keeps its initial value (Ca).
2023	ELSE
2024	! Eq. 24 of Yin et al. (2009)
2025	Obs = ( alpha(jv) * assimi(icinic,jl) ) / ( 0.047 * gbs(jv) ) + Oi
2026	! Eq. 23 of Yin et al. (2009)
2027	Cc(icinic,jl) = ( ( assimi(icinic,jl) + Rd(icinic,jl) ) * ( x2 * Obs + x3 ) + low_gamma_star(icinic) &
2028	* Obs * x1 ) &
2029	/ MAX(min_sechiba, x1 - ( assimi(icinic,jl) + Rd(icinic,jl) ))
2030	! Eq. 22 of Yin et al. (2009)
2031	leaf_ci(icinic,jv,jl) = ( Cc(icinic,jl) - ( b - assimi(icinic,jl) - Rm ) / gbs(jv) ) / a
2032	! Eq. 25 of Yin et al. (2009)
2033	! It should be Cs instead of Ca but it seems that
2034	! other equations in Appendix C make use of Ca
2035	gs(icinic,jl) = g0var(icinic) + ( assimi(icinic,jl) + Rd(icinic,jl) ) / &
2036	( Ca(icinic) - Cs_star ) * fvpd(icinic)
2037	ENDIF
2038	ENDDO
2039	ENDDO
2040	ENDIF
2041	ELSE
2042	!
2043	! @addtogroup Photosynthesis
2044	! @{
2045	!
2046	! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
2047	!! \latexonly
2048	!! \input{diffuco_trans_co2_2.4.3.tex}
2049	!! \endlatexonly
2050	! @}
2051	!
2052	!
2053	IF (nic .GT. 0) THEN
2054	DO inic=1,nic
2055	icinic=index_calc(inic)
2056
2057	A_1=9999.
2058
2059	! See eq. right after eq. 18 of Yin et al. (2009)
2060	DO limit_photo=1,2
2061	! Is Vc limiting the Assimilation
2062	IF ( limit_photo .EQ. 1 ) THEN
2063	x1 = vc2(icinic,jl)
2064	! It should be O not Oi (comment from Vuichard)
2065	x2 = KmC(icinic) * ( 1. + 2gamma_star(icinic)Sco(icinic) / KmO(icinic) )
2066	! Is J limiting the Assimilation
2067	ELSE
2068	x1 = JJ(icinic,jl)/4.
2069	x2 = 2. * gamma_star(icinic)
2070	ENDIF
2071
2072
2073	! See Appendix B of Yin et al. (2009)
2074	a = g0var(icinic) * ( x2 + gamma_star(icinic) ) + &
2075	( g0var(icinic) / gm(icinic) + fvpd(icinic) ) * ( x1 - Rd(icinic,jl) )
2076	b = Ca(icinic) * ( x1 - Rd(icinic,jl) ) - gamma_star(icinic) * x1 - Rd(icinic,jl) * x2
2077	c = Ca(icinic) + x2 + ( 1./gm(icinic) + 1./gb_co2(icinic) ) * ( x1 - Rd(icinic,jl) )
2078	d = x2 + gamma_star(icinic) + ( x1 - Rd(icinic,jl) ) / gm(icinic)
2079	m = 1./gm(icinic) + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * ( 1./gm(icinic) + 1./gb_co2(icinic) )
2080
2081	p = - ( d + (x1 - Rd(icinic,jl) ) / gm(icinic) + a * ( 1./gm(icinic) + 1./gb_co2(icinic) ) + &
2082	( g0var(icinic)/gm(icinic) + fvpd(icinic) ) * c ) / m
2083
2084	q = ( d * ( x1 - Rd(icinic,jl) ) + ac + ( g0var(icinic)/gm(icinic) + fvpd(icinic) ) b ) / m
2085	r = - a * b / m
2086
2087	! See Yin et al. (2009)
2088	QQ = ( (p*2._r_std) - 3._r_std q) / 9._r_std
2089	UU = ( 2._r_std* (p*3._r_std) - 9._r_std pq + 27._r_std r) /54._r_std
2090
2091	IF ( (QQ .GE. 0._r_std) .AND. (ABS(UU/(QQ**1.5_r_std) ) .LE. 1._r_std) ) THEN
2092	PSI = ACOS(UU/(QQ**1.5_r_std))
2093	A_1_tmp = -2._r_std * SQRT(QQ) * COS( PSI / 3._r_std ) - p / 3._r_std
2094	IF (( A_1_tmp .LT. A_1 )) THEN
2095	A_1 = A_1_tmp
2096	info_limitphoto(icinic,jl)=2.
2097	ELSE
2098	! In case, J is not limiting the assimilation
2099	! we have to re-initialise x1 and x2 values
2100	! in agreement with a Vc-limited assimilation
2101	x1 = vc2(icinic,jl)
2102	! It should be O not Oi (comment from Vuichard)
2103	x2 = KmC(icinic) * ( 1. + 2gamma_star(icinic)Sco(icinic) / KmO(icinic) )
2104	info_limitphoto(icinic,jl)=1.
2105	ENDIF
2106	ENDIF
2107	ENDDO
2108	IF ( (A_1 .EQ. 9999.) .OR. ( A_1 .LT. (-Rd(icinic,jl)) ) ) THEN
2109	IF ( printlev>=4 ) THEN
2110	WRITE(numout,*) 'We have a problem in diffuco_trans_co2 for A_1'
2111	WRITE(numout,*) 'no real positive solution found for pft:',jv
2112	WRITE(numout,*) 'temp_air:',temp_air(icinic)
2113	WRITE(numout,*) 'vpd:',vpd(icinic)
2114	END IF
2115	A_1 = -Rd(icinic,jl)
2116	ENDIF
2117	assimi(icinic,jl) = A_1
2118
2119	IF ( ABS( assimi(icinic,jl) + Rd(icinic,jl) ) .LT. min_sechiba ) THEN
2120	gs(icinic,jl) = g0var(icinic)
2121	ELSE
2122	! Eq. 18 of Yin et al. (2009)
2123	Cc(icinic,jl) = ( gamma_star(icinic) * x1 + ( assimi(icinic,jl) + Rd(icinic,jl) ) * x2 ) &
2124	/ MAX( min_sechiba, x1 - ( assimi(icinic,jl) + Rd(icinic,jl) ) )
2125	! Eq. 17 of Yin et al. (2009)
2126	leaf_ci(icinic,jv,jl) = Cc(icinic,jl) + assimi(icinic,jl) / gm(icinic)
2127	! See eq. right after eq. 15 of Yin et al. (2009)
2128	ci_star = gamma_star(icinic) - Rd(icinic,jl) / gm(icinic)
2129	!
2130	! Eq. 15 of Yin et al. (2009)
2131	gs(icinic,jl) = g0var(icinic) + ( assimi(icinic,jl) + Rd(icinic,jl) ) / ( leaf_ci(icinic,jv,jl) &
2132	- ci_star ) * fvpd(icinic)
2133	ENDIF
2134	ENDDO
2135	ENDIF
2136	ENDIF
2137	!
2138	IF (nic .GT. 0) THEN
2139	!
2140	DO inic=1,nic
2141	!
2142	! @addtogroup Photosynthesis
2143	! @{
2144	!
2145	!! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n
2146	!! \latexonly
2147	!! \input{diffuco_trans_co2_2.4.4.tex}
2148	!! \endlatexonly
2149	! @}
2150	!
2151	icinic=index_calc(inic)
2152	!
2153	! keep stomatal conductance of topmost level
2154	!
2155	IF ( jl .EQ. 1 ) THEN
2156	leaf_gs_top(icinic) = gs(icinic,jl)
2157	!
2158	ENDIF
2159	!
2160	! @addtogroup Photosynthesis
2161	! @{
2162	!
2163	!! 2.4.5 Integration at the canopy level\n
2164	!! \latexonly
2165	!! \input{diffuco_trans_co2_2.4.5.tex}
2166	!! \endlatexonly
2167	! @}
2168	! total assimilation and conductance
2169	assimtot(icinic) = assimtot(icinic) + &
2170	assimi(icinic,jl) * (laitab(jl+1)-laitab(jl))
2171	Rdtot(icinic) = Rdtot(icinic) + &
2172	Rd(icinic,jl) * (laitab(jl+1)-laitab(jl))
2173	gstot(icinic) = gstot(icinic) + &
2174	gs(icinic,jl) * (laitab(jl+1)-laitab(jl))
2175	!
2176	ilai(icinic) = jl
2177	!
2178	ENDDO
2179	!
2180	ENDIF
2181	ENDDO ! loop over LAI steps
2182
2183	IF(jv==testpft) THEN
2184	templeafci(:,:)=leaf_ci(:,testpft,:)
2185	CALL histwrite_p(hist_id, 'Cc', kjit, Cc, kjpindex*(nlai+1), indexlai)
2186	CALL histwrite_p(hist_id, 'Vc', kjit, Vc2, kjpindex*(nlai+1), indexlai)
2187	CALL histwrite_p(hist_id, 'Vj', kjit, JJ, kjpindex*(nlai+1), indexlai)
2188	CALL histwrite_p(hist_id, 'limitphoto', kjit, info_limitphoto, kjpindex*(nlai+1), indexlai)
2189	CALL histwrite_p(hist_id, 'gammastar', kjit, gamma_star, kjpindex,index)
2190	CALL histwrite_p(hist_id, 'Kmo', kjit, Kmo, kjpindex,index)
2191	CALL histwrite_p(hist_id, 'Kmc', kjit, Kmc, kjpindex,index)
2192	CALL histwrite_p(hist_id, 'gm', kjit, gm, kjpindex, index)
2193	CALL histwrite_p(hist_id, 'gs', kjit, gs, kjpindex*(nlai+1), indexlai)
2194	CALL histwrite_p(hist_id, 'leafci', kjit, templeafci, kjpindex*(nlai), indexlai)
2195	CALL histwrite_p(hist_id, 'assimi', kjit, assimi, kjpindex*(nlai+1), indexlai)
2196	CALL histwrite_p(hist_id, 'Rd', kjit, Rd, kjpindex*(nlai+1), indexlai)
2197	ENDIF
2198	!! Calculated intercellular CO2 over nlai needed for the chemistry module
2199	cim(:,jv)=0.
2200	laisum(:)=0
2201	DO jl=1,nlai
2202	WHERE (laitab(jl) .LE. lai(:,jv) )
2203	cim(:,jv)= cim(:,jv)+leaf_ci(:,jv,jl)*(laitab(jl+1)-laitab(jl))
2204	laisum(:)=laisum(:)+ (laitab(jl+1)-laitab(jl))
2205	ENDWHERE
2206	ENDDO
2207	WHERE (laisum(:)>0)
2208	cim(:,jv)= cim(:,jv)/laisum(:)
2209	ENDWHERE
2210
2211
2212	!
2213	!! 2.5 Calculate resistances
2214	!
2215	IF (nia .GT. 0) THEN
2216	!
2217	DO inia=1,nia
2218	!
2219	iainia=index_assi(inia)
2220
2221	!! Mean stomatal conductance for CO2 (mol m-2 s-1)
2222	gsmean(iainia,jv) = gstot(iainia)
2223	!
2224	! cimean is the "mean ci" calculated in such a way that assimilation
2225	! calculated in enerbil is equivalent to assimtot
2226	!
2227	IF ( ABS(gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) .GT. min_sechiba) THEN
2228	cimean(iainia,jv) = (fvpd(iainia)*(assimtot(iainia)+Rdtot(iainia))) /&
2229	(gsmean(iainia,jv)-g0var(iainia)*laisum(iainia)) + gamma_star(iainia)
2230	ELSE
2231	cimean(iainia,jv) = gamma_star(iainia)
2232	ENDIF
2233
2234	! conversion from umol m-2 (PFT) s-1 to gC m-2 (mesh area) tstep-1
2235	gpp(iainia,jv) = assimtot(iainia)12e-6veget_max(iainia,jv)*dt_sechiba
2236
2237	!
2238	! conversion from mol/m^2/s to m/s
2239	!
2240	! As in Pearcy, Schulze and Zimmermann
2241	! Measurement of transpiration and leaf conductance
2242	! Chapter 8 of Plant Physiological Ecology
2243	! Field methods and instrumentation, 1991
2244	! Editors:
2245	!
2246	! Robert W. Pearcy,
2247	! James R. Ehleringer,
2248	! Harold A. Mooney,
2249	! Philip W. Rundel
2250	!
2251	! ISBN: 978-0-412-40730-7 (Print) 978-94-010-9013-1 (Online)
2252
2253	gstot(iainia) = mol_to_m_1 (temp_air(iainia)/tp_00)&
2254	(pb_std/pb(iainia))gstot(iainia)ratio_H2O_to_CO2
2255	gstop(iainia) = mol_to_m_1 * (temp_air(iainia)/tp_00)*&
2256	(pb_std/pb(iainia))leaf_gs_top(iainia)ratio_H2O_to_CO2*&
2257	laitab(ilai(iainia)+1)
2258	!
2259	rveget(iainia,jv) = un/gstop(iainia)
2260
2261	!
2262	!
2263	! rstruct is the difference between rtot (=1./gstot) and rveget
2264	!
2265	! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
2266	!rstruct(iainia,jv) = un/gstot(iainia) - &
2267	! rveget(iainia,jv)
2268	rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
2269	rveget(iainia,jv), min_sechiba)
2270	!
2271	!
2272	!! wind is a global variable of the diffuco module.
2273	speed = MAX(min_wind, wind(iainia))
2274	!
2275	! beta for transpiration
2276	!
2277	! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
2278	!! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
2279	!! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
2280	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2281	! (un - zqsvegrap(iainia)) * &
2282	! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
2283	! rstruct(iainia,jv))))
2284	!! Global resistance of the canopy to evaporation
2285	cresist=(un / (un + speed * q_cdrag(iainia) * &
2286	veget(iainia,jv)/veget_max(iainia,jv) * &
2287	(rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
2288
2289	IF ( humrel(iainia,jv) >= min_sechiba ) THEN
2290	vbeta3(iainia,jv) = veget(iainia,jv) * &
2291	(un - zqsvegrap(iainia)) * cresist + &
2292	MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
2293	zqsvegrap(iainia) * cresist )
2294	ELSE
2295	! Because of a minimum conductance g0, vbeta3 cannot be zero even if humrel=0
2296	! in the above equation.
2297	! Here, we force transpiration to be zero when the soil cannot deliver it
2298	vbeta3(iainia,jv) = zero
2299	END IF
2300
2301	! vbeta3pot for computation of potential transpiration (needed for irrigation)
2302	vbeta3pot(iainia,jv) = MAX(zero, veget(iainia,jv) * cresist)
2303	!
2304	!
2305	ENDDO
2306	!
2307	ENDIF
2308	!
2309	END DO ! loop over vegetation types
2310
2311	IF (printlev>=3) WRITE (numout,*) ' diffuco_trans_co2 done '
2312
2313
2314	END SUBROUTINE diffuco_trans_co2
2315
2316
2317	!! ================================================================================================================================
2318	!! SUBROUTINE : diffuco_comb
2319	!!
2320	!>\BRIEF This routine combines the previous partial beta
2321	!! coefficients and calculates the total alpha and complete beta coefficients.
2322	!!
2323	!! DESCRIPTION : Those integrated coefficients are used to calculate (in enerbil.f90) the total evapotranspiration
2324	!! from the grid-cell. \n
2325	!!
2326	!! In the case that air is more humid than surface, dew deposition can occur (negative latent heat flux).
2327	!! In this instance, for temperature above zero, all of the beta coefficients are set to 0, except for
2328	!! interception (vbeta2) and bare soil (vbeta4 with zero soil resistance). The amount of water that is
2329	!! intercepted by leaves is calculated based on the value of LAI of the surface. In the case of freezing
2330	!! temperatures, water is added to the snow reservoir, and so vbeta4 and vbeta2 are set to 0, and the
2331	!! total vbeta is set to 1.\n
2332	!!
2333	!! \latexonly
2334	!! \input{diffucocomb1.tex}
2335	!! \endlatexonly
2336	!!
2337	!! The beta and alpha coefficients are initially set to 1.
2338	!! \latexonly
2339	!! \input{diffucocomb2.tex}
2340	!! \endlatexonly
2341	!!
2342	!! If snow is lower than the critical value:
2343	!! \latexonly
2344	!! \input{diffucocomb3.tex}
2345	!! \endlatexonly
2346	!! If in the presence of dew:
2347	!! \latexonly
2348	!! \input{diffucocomb4.tex}
2349	!! \endlatexonly
2350	!!
2351	!! Determine where the water goes (soil, vegetation, or snow)
2352	!! when air moisture exceeds saturation.
2353	!! \latexonly
2354	!! \input{diffucocomb5.tex}
2355	!! \endlatexonly
2356	!!
2357	!! If it is not freezing dew is put into the interception reservoir and onto the bare soil. If it is freezing,
2358	!! water is put into the snow reservoir.
2359	!! Now modify vbetas where necessary: for soil and snow
2360	!! \latexonly
2361	!! \input{diffucocomb6.tex}
2362	!! \endlatexonly
2363	!!
2364	!! and for vegetation
2365	!! \latexonly
2366	!! \input{diffucocomb7.tex}
2367	!! \endlatexonly
2368	!!
2369	!! Then compute part of dew that can be intercepted by leafs.
2370	!!
2371	!! There will be no transpiration when air moisture is too high, under any circumstance
2372	!! \latexonly
2373	!! \input{diffucocomb8.tex}
2374	!! \endlatexonly
2375	!!
2376	!! There will also be no interception loss on bare soil, under any circumstance.
2377	!! \latexonly
2378	!! \input{diffucocomb9.tex}
2379	!! \endlatexonly
2380	!!
2381	!! The flowchart details the 'decision tree' which underlies the module.
2382	!!
2383	!! RECENT CHANGE(S): None
2384	!!
2385	!! MAIN OUTPUT VARIABLE(S): vbeta1, vbeta4, humrel, vbeta2, vbeta3, vbeta
2386	!!
2387	!! REFERENCE(S) :
2388	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
2389	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
2390	!! Circulation Model. Journal of Climate, 6, pp.248-273
2391	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
2392	!! sur le cycle de l'eau, PhD Thesis, available from:
2393	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
2394	!!
2395	!! FLOWCHART :
2396	!! \latexonly
2397	!! \includegraphics[scale=0.25]{diffuco_comb_flowchart.png}
2398	!! \endlatexonly
2399	!! \n
2400	!_ ================================================================================================================================
2401
2402	SUBROUTINE diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
2403	& snow, veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, &
2404	& evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax)
2405
2406	! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006
2407
2408	!! 0. Variable and parameter declaration
2409
2410	!! 0.1 Input variables
2411
2412	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
2413	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
2414	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
2415	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Nortward Lowest level wind speed (m s^{-1})
2416	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
2417	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (hPa)
2418	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
2419	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
2420	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lower air temperature (K)
2421	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
2422	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (fraction)
2423	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
2424	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
2425	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil!! Total evaporating bare soil fraction
2426	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
2427
2428	!! 0.2 Output variables
2429
2430	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
2431
2432	!! 0.3 Modified variables
2433
2434	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation (-)
2435	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation (-)
2436	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
2437	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for interception loss (-)
2438	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration (-)
2439	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: evap_bare_lim !! limiting factor for bare soil evaporation
2440	!! when the 11-layer hydrology is used (-)
2441	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout):: evap_bare_lim_ns !! limiting factor for bare soil evaporation
2442	!! when the 11-layer hydrology is used (-)
2443	!! 0.4 Local variables
2444
2445	INTEGER(i_std) :: ji, jv
2446	REAL(r_std) :: zevtest, zsoil_moist, zrapp
2447	REAL(r_std), DIMENSION(kjpindex) :: qsatt
2448	LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow
2449	REAL(r_std) :: coeff_dew_veg
2450	REAL(r_std), DIMENSION(kjpindex) :: vegtot
2451	!_ ================================================================================================================================
2452
2453	!! 1 If we are in presence of dew
2454
2455	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
2456
2457
2458	!! 1.1 Determine where the water goes
2459	!! Determine where the water goes (soil, vegetation, or snow)
2460	!! when air moisture exceeds saturation.
2461	!! \latexonly
2462	!! \input{diffucocomb5.tex}
2463	!! \endlatexonly
2464	toveg(:) = .FALSE.
2465	tosnow(:) = .FALSE.
2466	DO ji = 1, kjpindex
2467	IF ( qsatt(ji) .LT. qair(ji) ) THEN
2468	IF (temp_air(ji) .GT. tp_00) THEN
2469	!! If it is not freezing dew is put into the
2470	!! interception reservoir and onto the bare soil.
2471	toveg(ji) = .TRUE.
2472	ELSE
2473	!! If it is freezing water is put into the
2474	!! snow reservoir.
2475	tosnow(ji) = .TRUE.
2476	ENDIF
2477	ENDIF
2478	END DO
2479
2480	!! 1.2 Now modify vbetas where necessary.
2481
2482	!! 1.2.1 Soil and snow
2483	!! \latexonly
2484	!! \input{diffucocomb6.tex}
2485	!! \endlatexonly
2486
2487	! We need to keep consistency between evap_bare_lim, evap_bare_lim_ns and vbeta4 (thus vevapnu)
2488	! or we have a water conservation issue in hydrol_split_soil
2489
2490	DO ji = 1, kjpindex
2491
2492	IF ( toveg(ji) ) THEN
2493
2494	vbeta1(ji) = zero
2495	vegtot(ji) = SUM(veget_max(ji,:))
2496
2497	IF ( (tot_bare_soil(ji) .GT. min_sechiba) .AND. (vegtot(ji).GT. min_sechiba) ) THEN
2498
2499	vbeta4(ji) = tot_bare_soil(ji)
2500
2501	! We now have to redefine evap_bare_lim & evap_bare_lim_ns
2502	IF (evap_bare_lim(ji) .GT. min_sechiba) THEN
2503	evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji)
2504	ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles
2505	evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji)
2506	ENDIF
2507
2508	evap_bare_lim(ji) = vbeta4(ji)
2509	! consistent with evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)soiltile(ji,:)vegtot(ji))
2510	! as SUM(soiltile(ji,:)) = 1
2511
2512	ELSE
2513	vbeta4(ji) = zero
2514	evap_bare_lim_ns(ji,:) = zero
2515	evap_bare_lim(ji) = zero
2516	ENDIF
2517	ENDIF
2518
2519	IF ( tosnow(ji) ) THEN
2520	vbeta1(ji) = un
2521	vbeta4(ji) = zero
2522	evap_bare_lim_ns(ji,:) = zero
2523	evap_bare_lim(ji) = zero
2524	ENDIF
2525
2526	ENDDO
2527
2528	!! 1.2.2 Vegetation and interception loss
2529	!! \latexonly
2530	!! \input{diffucocomb7.tex}
2531	!! \endlatexonly
2532	DO jv = 1, nvm
2533
2534	DO ji = 1, kjpindex
2535
2536	IF ( toveg(ji) ) THEN
2537	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
2538
2539	! Compute part of dew that can be intercepted by leafs.
2540	IF ( lai(ji,jv) .GT. min_sechiba) THEN
2541	IF (lai(ji,jv) .GT. 1.5) THEN
2542	coeff_dew_veg= &
2543	& dew_veg_poly_coeff(6)lai(ji,jv)*5 &
2544	& - dew_veg_poly_coeff(5)lai(ji,jv)*4 &
2545	& + dew_veg_poly_coeff(4)lai(ji,jv)*3 &
2546	& - dew_veg_poly_coeff(3)lai(ji,jv)*2 &
2547	& + dew_veg_poly_coeff(2)*lai(ji,jv) &
2548	& + dew_veg_poly_coeff(1)
2549	ELSE
2550	coeff_dew_veg=un
2551	ENDIF
2552	ELSE
2553	coeff_dew_veg=zero
2554	ENDIF
2555	IF (jv .EQ. 1) THEN
2556	! This line may not work with CWRR when frac_bare is distributed among three soiltiles
2557	! Fortunately, qsintmax(ji,1)=0 (LAI=0 in PFT1) so we never pass here
2558	vbeta2(ji,jv) = coeff_dew_veg*tot_bare_soil(ji)
2559	ELSE
2560	vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv)
2561	ENDIF
2562	ELSE
2563	vbeta2(ji,jv) = zero ! if qsintmax=0, vbeta2=0
2564	ENDIF
2565	ENDIF
2566	IF ( tosnow(ji) ) vbeta2(ji,jv) = zero
2567
2568	ENDDO
2569
2570	ENDDO
2571
2572	!! 1.2.3 Vegetation and transpiration
2573	!! There will be no transpiration when air moisture is too high, under any circumstance
2574	!! \latexonly
2575	!! \input{diffucocomb8.tex}
2576	!! \endlatexonly
2577	DO jv = 1, nvm
2578	DO ji = 1, kjpindex
2579	IF ( qsatt(ji) .LT. qair(ji) ) THEN
2580	vbeta3(ji,jv) = zero
2581	humrel(ji,jv) = zero
2582	ENDIF
2583	ENDDO
2584	ENDDO
2585
2586
2587	!! 1.2.4 Overrules 1.2.2
2588	!! There will also be no interception loss on bare soil, under any circumstance.
2589	!! \latexonly
2590	!! \input{diffucocomb9.tex}
2591	!! \endlatexonly
2592	DO ji = 1, kjpindex
2593	IF ( qsatt(ji) .LT. qair(ji) ) THEN
2594	vbeta2(ji,1) = zero
2595	ENDIF
2596	ENDDO
2597
2598	!! 2 Now calculate vbeta in all cases (the equality needs to hold for enerbil to be consistent)
2599
2600	DO ji = 1, kjpindex
2601	vbeta(ji) = vbeta4(ji) + SUM(vbeta2(ji,:)) + SUM(vbeta3(ji,:))
2602
2603	IF (vbeta(ji) .LT. min_sechiba) THEN
2604	vbeta(ji) = zero
2605	vbeta4(ji) = zero
2606	vbeta2(ji,:)= zero
2607	vbeta3(ji,:)= zero
2608	evap_bare_lim_ns(ji,:) = zero
2609	evap_bare_lim(ji) = zero
2610	END IF
2611	ENDDO
2612
2613	CALL xios_orchidee_send_field("evap_bare_lim",evap_bare_lim)
2614	CALL xios_orchidee_send_field("evap_bare_lim_ns",evap_bare_lim_ns)
2615
2616	IF (printlev>=3) WRITE (numout,*) ' diffuco_comb done '
2617
2618	END SUBROUTINE diffuco_comb
2619
2620
2621	!! ================================================================================================================================
2622	!! SUBROUTINE : diffuco_raerod
2623	!!
2624	!>\BRIEF Computes the aerodynamic resistance, for cases in which the
2625	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ and when the flag
2626	!! 'ldq_cdrag_from_gcm' is set to TRUE
2627	!!
2628	!! DESCRIPTION : Simply computes the aerodynamic resistance, for cases in which the
2629	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient
2630	!! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine
2631	!! diffuco_aero is called instead of this one.
2632	!!
2633	!! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed:
2634	!! \latexonly
2635	!! \input{diffucoaerod1.tex}
2636	!! \endlatexonly
2637	!!
2638	!! next calculate ::raero
2639	!! \latexonly
2640	!! \input{diffucoaerod2.tex}
2641	!! \endlatexonly
2642	!!
2643	!! RECENT CHANGE(S): None
2644	!!
2645	!! MAIN OUTPUT VARIABLE(S): ::raero
2646	!!
2647	!! REFERENCE(S) :
2648	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
2649	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
2650	!! Circulation Model. Journal of Climate, 6, pp.248-273
2651	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influence de l'irrigation
2652	!! sur le cycle de l'eau, PhD Thesis, available from:
2653	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
2654	!!
2655	!! FLOWCHART : None
2656	!! \n
2657	!_ ================================================================================================================================
2658
2659	SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
2660
2661	IMPLICIT NONE
2662
2663	!! 0. Variable and parameter declaration
2664
2665	!! 0.1 Input variables
2666
2667	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
2668	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind velocity (m s^{-1})
2669	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind velocity (m s^{-1})
2670	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag coefficient (-)
2671
2672	!! 0.2 Output variables
2673
2674	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance (s m^{-1})
2675
2676	!! 0.3 Modified variables
2677
2678	!! 0.4 Local variables
2679
2680	INTEGER(i_std) :: ji !! (-)
2681	REAL(r_std) :: speed !! (m s^{-1})
2682	!_ ================================================================================================================================
2683
2684	!! 1. Simple calculation of the aerodynamic resistance, for diganostic purposes.
2685
2686	DO ji=1,kjpindex
2687
2688	!! \latexonly
2689	!! \input{diffucoaerod1.tex}
2690	!! \endlatexonly
2691	speed = MAX(min_wind, wind(ji))
2692
2693	!! \latexonly
2694	!! \input{diffucoaerod2.tex}
2695	!! \endlatexonly
2696	raero(ji) = un / (q_cdrag(ji)*speed)
2697
2698	ENDDO
2699
2700	END SUBROUTINE diffuco_raerod
2701
2702
2703	FUNCTION Arrhenius (kjpindex,temp,ref_temp,energy_act) RESULT ( val_arrhenius )
2704	!! 0.1 Input variables
2705
2706	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
2707	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
2708	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
2709	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
2710
2711	!! 0.2 Result
2712
2713	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
2714	!! a Arrhenius function (-)
2715
2716	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:))))
2717	END FUNCTION Arrhenius
2718
2719	FUNCTION Arrhenius_modified_1d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
2720	!! 0.1 Input variables
2721
2722	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
2723	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
2724	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
2725	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
2726	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
2727	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
2728
2729	!! 0.2 Result
2730
2731	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
2732	!! a Arrhenius function (-)
2733
2734	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
2735	* (1. + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * RR ))) &
2736	/ (1. + EXP( (temp(:) * entropy(:) - energy_deact) / ( RR*temp(:))))
2737
2738	END FUNCTION Arrhenius_modified_1d
2739
2740	FUNCTION Arrhenius_modified_0d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
2741	!! 0.1 Input variables
2742
2743	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
2744	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
2745	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
2746	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
2747	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
2748	REAL(r_std),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
2749
2750	!! 0.2 Result
2751
2752	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
2753	!! a Arrhenius function (-)
2754
2755	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
2756	* (1. + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * RR ))) &
2757	/ (1. + EXP( (temp(:) * entropy - energy_deact) / ( RR*temp(:))))
2758
2759	END FUNCTION Arrhenius_modified_0d
2760
2761
2762	END MODULE diffuco

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

Download in other formats: