Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_sechiba/diffuco.f90 @ 7346

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

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

Download in other formats: