Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_sechiba/diffuco.f90

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

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

Download in other formats: