Context Navigation

diffuco.f90 @ 7709

Last change on this file since 7709 was 7709, checked in by josefine.ghattas, 2 years ago

Integrated new irrigation scheme developed by Pedro Arboleda. See ticket #857
This corresponds to revsion 7708 of version pedro.arboleda/ORCHIDEE. Following differences were made but were not made on the pedro.arboleda/ORCHIDEE :

argumet place in call to routing_wrapper_intialize changed order
lines with only change in space were not taken
some indentation changed
set irrigation output as enalbled false if not do_irrigation

Property svn:keywords set to Revision Date HeadURL Date Author Revision

File size: 136.1 KB

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

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

Context Navigation

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

Download in other formats: