Context Navigation

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

Last change on this file since 7251 was 7251, checked in by agnes.ducharne, 3 years ago
Same correction as in r6465 in the trunk: Correct the adjustments to evap_bare_lim_ns in diffuco_bare.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 135.0 KB

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

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

Download in other formats: