Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/diffuco.f90 @ 5491

Last change on this file since 5491 was 5488, checked in by chunjing.qiu, 6 years ago
C balance checked
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 153.9 KB

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

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

Download in other formats: