Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diffuco.f90 @ 7541

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/diffuco.f90 @ 7541

Download in other formats: