Context Navigation

diffuco.f90 @ 854

Last change on this file since 854 was 764, checked in by james.ryder, 12 years ago
Revision for diffuco.f90 JR 190312
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 122.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	!! - diffuco_inca 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_veg
56	USE sechiba_io
57	USE ioipsl
58	USE constantes_co2
59	USE parallel
60	! USE WRITE_FIELD_p
61
62	IMPLICIT NONE
63
64	! public routines :
65	! diffuco_main only
66	PRIVATE
67	PUBLIC :: diffuco_main,diffuco_clear
68
69	!
70	! variables used inside diffuco module : declaration and initialisation
71	!
72
73	! variables common to diffuco module
74	INTEGER(i_std), PARAMETER :: nlai = 20
75	LOGICAL, SAVE :: l_first_diffuco = .TRUE. !! Initialisation has to be done one time
76	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (-)
77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: leaf_ci !! intercellular CO2 concentration (ppm)
78	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rstruct !! architectural resistance (s m^{-1})
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: raero !! Aerodynamic resistance (s m^{-1})
80	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsatt !! Surface saturated humidity (kg kg^{-1})
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind module (m s^{-1})
82	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rveg_pft !! Nathalie le 28 mars 2006 - sur proposition de Fred Hourdin, ajout
83	!! d'un potentiometre pour regler la resistance de la vegetation
84
85	CONTAINS
86
87
88	!! ================================================================================================================================
89	!! SUBROUTINE : diffuco_main
90	!!
91	!>\BRIEF The root subroutine for the module, which calls all other required
92	!! subroutines.
93	!!
94	!! DESCRIPTION :
95
96	!! This is the root subroutine for the module. Following
97	!! initialisation (if required) and preparation of the restart file, the module first of all calculates
98	!! the surface drag coefficient (via a call to diffuco_aero), using available parameters to determine
99	!! the stability of air in the surface layer by calculating the Richardson Nubmber. If a value for the
100	!! surface drag coefficient is passed down from the atmospheric model and and if the flag 'ldq_cdrag_from_gcm'
101	!! is set to TRUE, then the subroutine diffuco_aerod is called instead. This calculates the aerodynamic coefficient. \n
102	!!
103	!! Following this, an estimation of the saturated humidity at the surface is made (via a call
104	!! to qsatcalc in the module qsat_moisture). Following this the beta coefficients for sublimation (via
105	!! diffuco_snow), interception (diffuco_inter), bare soil (diffuco_bare), and transpiration (via
106	!! diffuco_trans_co2 if co2 is considered, diffuco_trans otherwise) are calculated in sequence. Finally
107	!! the alpha and beta coefficients are combined (diffuco_comb). \n
108	!!
109	!! The surface drag coefficient is calculated for use within the module enerbil. It is required to to
110	!! calculate the aerodynamic coefficient for every flux. \n
111	!!
112	!! The various beta coefficients are used within the module enerbil for modifying the rate of evaporation,
113	!! as appropriate for the surface. As explained in Chapter 2 of Guimberteau (2010), that module (enerbil)
114	!! calculates the rate of evaporation essentially according to the expression $E = /beta E_{pot}$, where
115	!! E is the total evaporation and $E_{pot}$ the evaporation potential. If $\beta = 1$, there would be
116	!! essentially no resistance to evaporation, whereas should $\beta = 0$, there would be no evaporation and
117	!! the surface layer would be subject to some very stong hydrological stress. \n
118	!!
119	!! ATTENTION ::valpha [DISPOSABLE] is not used anymore ! Its value is always 1.
120	!!
121	!! RECENT CHANGE(S): None
122	!!
123	!! MAIN OUTPUT VARIABLE(S): humrel, q_cdrag, vbeta, valpha, vbeta1, vbeta4,
124	!! vbetaco2, vbeta2, vbeta3, rveget, cimean
125	!!
126	!! REFERENCE(S) :
127	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
128	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
129	!! Circulation Model. Journal of Climate, 6, pp.248-273.
130	!! - de Rosnay, P, 1999. ReprÃ©sentation des interactions sol-plante-atmosphÃšre dans le modÃšle de circulation gÃ©nÃ©rale
131	!! du LMD, 1999. PhD Thesis, UniversitÃ© Paris 6, available (25/01/12):
132	!! http://www.ecmwf.int/staff/patricia_de_rosnay/publications.html#8
133	!! - Ducharne, A, 1997. Le cycle de l'eau: modÃ©lisation de l'hydrologie continentale, Ã©tude de ses interactions avec
134	!! le climat, PhD Thesis, UniversitÃ© Paris 6
135	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
136	!! sur le cycle de l'eau, PhD Thesis, available (25/01/12):
137	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
138	!! - LathiÃšre, J, 2005. Evolution des Ã©missions de composÃ©s organiques et azotÃ©s par la biosphÃšre continentale dans le
139	!! modÃšle LMDz-INCA-ORCHIDEE, UniversitÃ© Paris 6
140	!!
141	!! FLOWCHART :
142	!! \latexonly
143	!! \includegraphics[scale=0.5]{diffuco_main_flowchart.png}
144	!! \endlatexonly
145	!! \n
146	!_ ================================================================================================================================
147
148
149	! Main routine for diffuco module.
150	! - called only one time for initialisation
151	! - called every time step for calculations (also the first time step)
152	! - called one more time at last time step for writing _restart_ file
153	!
154	! The following processes are calculated:
155	! - call diffuco_aero for aerodynamic transfer coeficient
156	! - call diffuco_snow for partial beta coefficient: sublimation
157	! - call diffuco_inter for partial beta coefficient: interception for each type of vegetation
158	! - call diffuco_bare for partial beta coefficient: bare soil
159	! - call diffuco_trans for partial beta coefficient: transpiration for each type of vegetation, using Jarvis formula
160	! - call diffuco_trans_co2 for partial beta coefficient: transpiration for each type of vegetation, using Farqhar's formula
161	! - call diffuco_comb for alpha and beta coefficient
162	! - call diffuco_inca for alpha and beta coefficients for biogenic emissions
163	SUBROUTINE diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
164	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
165	! & zlev, z0, roughheight, temp_sol, temp_air, rau, q_cdrag, qsurf, qair, pb, &
166	& zlev, z0, roughheight, temp_sol, temp_air, rau, q_cdrag, qsurf, qair, q2m, t2m, pb, &
167	& rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
168	& swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
169	& vbeta , valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
170
171	!! 0. Variable and parameter declaration
172
173	!! 0.1 Input variables
174
175	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (-)
176	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
177	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier (-)
178	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier (-)
179	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier (-)
180	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
181	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart file to read (-)
182	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart file to write (-)
183	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
184	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map (-)
185	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
186	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
187	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer (m)
188	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0 !! Surface roughness Length (m)
189	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness (m)
190	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
191	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
192	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
193	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! Near surface air specific humidity (kg kg^{-1})
194	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level air specific humidity (kg kg^{-1})
195	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m air specific humidity (kg kg^{-1})
196	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature (K)
197	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
198	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure (Pa)
199	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rsol !! Bare soil evaporation resistance (s m^{-1})
200	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! Limit to the bare soil evaporation when the
201	!! 11-layer hydrology is used (-)
202	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (mm day^{-1})
203	!! NdN This variable does not seem to be used at
204	!! all in diffuco!
205	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,... (-)
206	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,... (kg m^{-2})
207	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... (-)
208	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux (W m^{-2})
209	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux (W m^{-2})
210	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
211	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (-)
212	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
213	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
214	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
215	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
216	!! (kg m^{-2})
217	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! min+max+opt temps, vcmax, vjmax
218	!! for photosynthesis (K ??)
219	!! 0.2 Output variables
220
221	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
222	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: valpha !! Total alpha coefficient (-)
223	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (-)
224	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
225	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE: Beta for CO2 assimilation (-)
226	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
227	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration (-)
228	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Stomatal resistance for the whole canopy (s m^{-1})
229	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE: mean intercellular ci (\mumole m^{-2} s^{-1})
230
231	!! 0.3 Modified variables
232
233	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
234	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Product of drag coefficient and wind speed (m s^{-1})
235
236	!! 0.4 Local variables
237
238	REAL(r_std),DIMENSION (kjpindex,nvm) :: vbeta23 !! Beta for fraction of wetted foliage that will
239	!! transpire once intercepted water has evaporated (-)
240	INTEGER(i_std) :: ilai
241	CHARACTER(LEN=4) :: laistring
242	!_ ================================================================================================================================
243
244	!! 1. Perform initialisation, if required
245
246	IF (l_first_diffuco) THEN
247
248	!Config Key = CDRAG_FROM_GCM
249	!Config Desc = Keep cdrag coefficient from gcm.
250	!Config Def = TRUE if q_cdrag on initialization is non zero
251	!Config Help = Set to .TRUE. if you want q_cdrag coming from GCM.
252	!Congig Keep cdrag coefficient from gcm for latent and sensible heat fluxes.
253	IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN
254	ldq_cdrag_from_gcm = .FALSE.
255	ELSE
256	ldq_cdrag_from_gcm = .TRUE.
257	ENDIF
258
259	!?? q_cdrag is always 0 on initialization ??
260	CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm)
261
262	WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm
263
264	IF (long_print) WRITE (numout,*) ' call diffuco_init '
265
266	! If cdrag is
267	CALL diffuco_init(kjit, ldrestart_read, kjpindex, index, rest_id, q_cdrag)
268
269	RETURN
270
271	ENDIF
272
273	!! 2. Prepare the restart file for the next simulation
274
275	IF (ldrestart_write) THEN
276
277	IF (long_print) WRITE (numout,*) ' we have to complete restart file with DIFFUCO variables '
278
279	var_name= 'rstruct'
280	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g)
281
282	var_name= 'raero'
283	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, raero, 'scatter', nbp_glo, index_g)
284
285	var_name= 'qsatt'
286	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, qsatt, 'scatter', nbp_glo, index_g)
287
288	! the following variable is written only if CO2 was calculated
289	IF ( control%ok_co2 ) THEN
290
291	DO ilai = 1, nlai
292
293	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
294	write(laistring,'(i4)') ilai
295	laistring=ADJUSTL(laistring)
296	var_name='leaf_ci_'//laistring(1:LEN_TRIM(laistring))
297
298	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, leaf_ci(:,:,ilai), 'scatter', nbp_glo, index_g)
299
300	ENDDO
301
302	ENDIF
303
304	IF (.NOT.ldq_cdrag_from_gcm) THEN
305	var_name= 'cdrag'
306	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, q_cdrag, 'scatter', nbp_glo, index_g)
307	ENDIF
308
309	RETURN
310
311	END IF
312
313	wind(:) = SQRT (u(:)u(:) + v(:)v(:))
314
315
316	!! 3. Calculate the different coefficients
317
318	IF (.NOT.ldq_cdrag_from_gcm) THEN
319
320	! Case 3a)
321	CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0, roughheight, veget_max, temp_sol, temp_air, &
322	& qsurf, qair, q_cdrag)
323
324	ENDIF
325
326	! Case 3b)
327	CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
328
329	!! 4. Make an estimation of the saturated humidity at the surface
330
331	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
332
333	!! 5. Calculate the beta coefficient for sublimation
334
335	CALL diffuco_snow (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, &
336	& snow, frac_nobio, totfrac_nobio, snow_nobio, vbeta1)
337
338	!! 6. Calculate the beta coefficient for interception
339
340	! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23
341	!CALL diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
342	! & qsintveg, qsintmax, rstruct, vbeta2)
343	CALL diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
344	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
345
346	!! 7. Calculate the beta coefficient for bare soil
347
348	CALL diffuco_bare (kjpindex, dtradia, u, v, q_cdrag, rsol, evap_bare_lim, evapot, humrel, veget, vbeta4)
349
350	!! 8. Calculate the beta coefficient for transpiration
351
352	IF ( control%ok_co2 ) THEN
353
354	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
355	! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23
356	!CALL diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
357	! assim_param, ccanopy, &
358	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
359
360	! case 8a)
361	CALL diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, &
362	assim_param, ccanopy, &
363	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
364
365	ELSE
366
367	! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23
368	!CALL diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
369	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
370
371	! case 8b)
372	CALL diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
373	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
374
375	ENDIF
376
377	!! 9. Combine the alpha and beta coefficients
378
379	! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006
380	CALL diffuco_comb (kjpindex, dtradia, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, &
381	& veget, lai, vbeta1, vbeta2, vbeta3, vbeta4, valpha, vbeta, qsintmax)
382
383	IF ( .NOT. almaoutput ) THEN
384	CALL histwrite(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
385	CALL histwrite(hist_id, 'raero', kjit, raero, kjpindex, index)
386	! Ajouts Nathalie - novembre 2006
387	CALL histwrite(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
388	CALL histwrite(hist_id, 'Wind', kjit, wind, kjpindex, index)
389	! Fin ajouts Nathalie
390	CALL histwrite(hist_id, 'qsatt', kjit, qsatt, kjpindex, index)
391
392	IF ( hist2_id > 0 ) THEN
393	CALL histwrite(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
394	CALL histwrite(hist2_id, 'raero', kjit, raero, kjpindex, index)
395	CALL histwrite(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
396	CALL histwrite(hist2_id, 'Wind', kjit, wind, kjpindex, index)
397	CALL histwrite(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index)
398	ENDIF
399	ELSE
400
401	ENDIF
402
403	IF (long_print) WRITE (numout,*) ' diffuco_main done '
404
405	END SUBROUTINE diffuco_main
406
407
408	!! ================================================================================================================================
409	!! SUBROUTINE : diffuco_init
410	!!
411	!>\BRIEF Dynamic allocation algorithm for local arrays
412	!!
413	!! DESCRIPTION : Housekeeping module to allocate local arrays
414	!!
415	!! RECENT CHANGE(S): None
416	!!
417	!! MAIN OUTPUT VARIABLE(S) : q_cdrag
418	!!
419	!! REFERENCE(S) : None
420	!!
421	!! FLOWCHART : None
422	!! \n
423	!_ ================================================================================================================================
424
425	SUBROUTINE diffuco_init(kjit, ldrestart_read, kjpindex, index, rest_id, q_cdrag)
426
427	!! 0. Variable and parameter declaration
428
429	!! 0.1 Input variables
430
431	INTEGER(i_std), INTENT (in) :: kjit !! Time step number (-)
432	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read (-)
433	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (-)
434	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map (-)
435	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier (-)
436	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Product of Surface drag and wind speed (m s^{-1})
437
438	!! 0.2 Output variables
439
440	!! 0.3 Modified variables
441
442	!! 0.4 Local variables
443
444	INTEGER(i_std) :: ier, jv
445	INTEGER(i_std) :: ilai
446	CHARACTER(LEN=4) :: laistring
447	REAL(r_std),DIMENSION (kjpindex) :: temp
448	!_ ================================================================================================================================
449
450	!! 1. Initialisation
451
452	IF (l_first_diffuco) THEN
453	l_first_diffuco=.FALSE.
454	ELSE
455	WRITE (numout,*) ' l_first_diffuco false . we stop '
456	STOP 'diffuco_init'
457	ENDIF
458
459	! allocate only if CO2 is calculated
460	IF ( control%ok_co2 ) THEN
461
462	ALLOCATE (leaf_ci(kjpindex,nvm,nlai),stat=ier)
463	IF (ier.NE.0) THEN
464	WRITE (numout,) ' error in leaf_ci allocation. We stop. We need kjpindexnvm*nlai words = ',&
465	kjpindexnvmnlai
466	STOP 'diffuco_init'
467	END IF
468
469	ENDIF
470
471	ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
472	IF (ier.NE.0) THEN
473	WRITE (numout,*) ' error in rstruct allocation. We stop. We need kjpindex x nvm words = ' ,kjpindex,' x ' ,nvm,&
474	& ' = ',kjpindex*nvm
475	STOP 'diffuco_init'
476	END IF
477	ALLOCATE (raero(kjpindex),stat=ier)
478	IF (ier.NE.0) THEN
479	WRITE (numout,*) ' error in raero allocation. We stop. We need kjpindex x nvm words = ', kjpindex
480	STOP 'diffuco_init'
481	END IF
482
483	ALLOCATE (qsatt(kjpindex),stat=ier)
484	IF (ier.NE.0) THEN
485	WRITE (numout,*) ' error in qsatt allocation. We stop. We need kjpindex x nvm words = ', kjpindex
486	STOP 'diffuco_init'
487	END IF
488
489	ALLOCATE (wind(kjpindex),stat=ier)
490	IF (ier.NE.0) THEN
491	WRITE (numout,*) ' error in wind allocation. We stop. We need kjpindex x nvm words = ', kjpindex
492	STOP 'diffuco_init'
493	END IF
494
495	IF (ldrestart_read) THEN
496
497	IF (long_print) WRITE (numout,*) ' we have to read a restart file for DIFFUCO variables'
498
499	var_name='rstruct'
500	CALL ioconf_setatt('UNITS', 's/m')
501	CALL ioconf_setatt('LONG_NAME','Structural resistance')
502	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g)
503
504	IF ( MINVAL(rstruct) .EQ. MAXVAL(rstruct) .AND. MAXVAL(rstruct) .EQ. val_exp ) THEN
505
506	DO jv = 1, nvm
507	rstruct(:,jv) = rstruct_const(jv)
508	ENDDO
509
510	ENDIF
511
512	var_name='raero' ;
513	CALL ioconf_setatt('UNITS', 's/m')
514	CALL ioconf_setatt('LONG_NAME','Aerodynamic resistance')
515	IF ( ok_var(var_name) ) THEN
516	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
517	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
518	raero(:) = temp(:)
519	ENDIF
520	ENDIF
521
522	var_name='qsatt' ;
523	CALL ioconf_setatt('UNITS', 'g/g')
524	CALL ioconf_setatt('LONG_NAME','Surface saturated humidity')
525	IF ( ok_var(var_name) ) THEN
526	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
527	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
528	qsatt(:) = temp(:)
529	ENDIF
530	ENDIF
531
532	! the following variable is read only if CO2 is calculated
533	IF ( control%ok_co2 ) THEN
534
535	CALL ioconf_setatt('UNITS', 'ppm')
536	CALL ioconf_setatt('LONG_NAME','Leaf CO2')
537
538	DO ilai = 1, nlai
539
540	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
541	write(laistring,'(i4)') ilai
542	laistring=ADJUSTL(laistring)
543	var_name='leaf_ci_'//laistring(1:LEN_TRIM(laistring))
544
545	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE.,leaf_ci(:,:,ilai), "gather", nbp_glo, index_g)
546
547	ENDDO
548
549	!Config Key = DIFFUCO_LEAFCI
550	!Config Desc = Initial leaf CO2 level if not found in restart
551	!Config Def = 233.
552	!Config Help = The initial value of leaf_ci if its value is not found
553	!Config in the restart file. This should only be used if the model is
554	!Config started without a restart file.
555	CALL setvar_p (leaf_ci, val_exp,'DIFFUCO_LEAFCI', 233._r_std)
556	ENDIF
557
558	IF (.NOT.ldq_cdrag_from_gcm) THEN
559	var_name= 'cdrag'
560	CALL ioconf_setatt('LONG_NAME','Drag coefficient for LE and SH')
561	CALL ioconf_setatt('UNITS', '-')
562	IF ( ok_var(var_name) ) THEN
563	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
564	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
565	q_cdrag(:) = temp(:)
566	ENDIF
567	ENDIF
568
569	ENDIF
570
571	ENDIF
572
573	!
574	! Ajouts Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin
575	!
576	ALLOCATE (rveg_pft(nvm),stat=ier)
577	IF (ier.NE.0) THEN
578	WRITE (numout,*) &
579	' error in rveg_pft allocation. We stop. We need nvm words = ', nvm
580	STOP 'diffuco_init'
581	END IF
582	!
583	!Config Key = RVEG_PFT
584	!Config Desc = Artificial parameter to increase or decrease canopy resistance.
585	!Config Def = 1.
586	!Config Help = This parameter is set by PFT.
587	rveg_pft(:) = un
588	CALL getin_p('RVEG_PFT', rveg_pft)
589
590	WRITE(numout,*) 'DANS DIFFUCO_INIT , RVEG_PFT=',rveg_pft
591
592	IF (long_print) WRITE (numout,*) ' diffuco_init done '
593
594	END SUBROUTINE diffuco_init
595
596
597	!! ================================================================================================================================
598	!! SUBROUTINE : diffuco_clear
599	!!
600	!>\BRIEF Housekeeping module to deallocate the variables
601	!! leaf_ci, rstruct and raero
602	!!
603	!! DESCRIPTION : Housekeeping module to deallocate the variables
604	!! leaf_ci, rstruct and raero
605	!!
606	!! RECENT CHANGE(S) : None
607	!!
608	!! MAIN OUTPUT VARIABLE(S) : None
609	!!
610	!! REFERENCE(S) : None
611	!!
612	!! FLOWCHART : None
613	!! \n
614	!_ ================================================================================================================================
615
616	SUBROUTINE diffuco_clear()
617
618	l_first_diffuco=.TRUE.
619
620	IF (ALLOCATED (leaf_ci)) DEALLOCATE (leaf_ci)
621	IF (ALLOCATED (rstruct)) DEALLOCATE (rstruct)
622	IF (ALLOCATED (raero)) DEALLOCATE (raero)
623	IF (ALLOCATED (rveg_pft)) DEALLOCATE (rveg_pft)
624
625	END SUBROUTINE diffuco_clear
626
627
628	!! ================================================================================================================================
629	!! SUBROUTINE : diffuco_aero
630	!!
631	!>\BRIEF This module first calculates the surface drag
632	!! coefficient, for cases in which the surface drag coefficient is NOT provided by the coupled
633	!! atmospheric model LMDZ or when the flag ldq_cdrag_from_gcm is set to FALSE
634	!!
635	!! DESCRIPTION : Computes the surface drag coefficient, for cases
636	!! in which it is NOT provided by the coupled atmospheric model LMDZ. The module first uses the
637	!! meteorolgical input to calculate the Richardson Number, which is an indicator of atmospheric
638	!! stability in the surface layer. The formulation used to find this surface drag coefficient is
639	!! dependent on the stability determined. \n
640	!!
641	!! Designation of wind speed
642	!! \latexonly
643	!! \input{diffucoaero1.tex}
644	!! \endlatexonly
645	!!
646	!! Calculation of geopotential. This is the definition of Geopotential height (e.g. Jacobson
647	!! eqn.4.47, 2005). (required for calculation of the Richardson Number)
648	!! \latexonly
649	!! \input{diffucoaero2.tex}
650	!! \endlatexonly
651	!!
652	!! \latexonly
653	!! \input{diffucoaero3.tex}
654	!! \endlatexonly
655	!!
656	!! Calculation of the virtual air temperature at the surface (required for calculation
657	!! of the Richardson Number)
658	!! \latexonly
659	!! \input{diffucoaero4.tex}
660	!! \endlatexonly
661	!!
662	!! Calculation of the virtual surface temperature (required for calculation of th
663	!! Richardson Number)
664	!! \latexonly
665	!! \input{diffucoaero5.tex}
666	!! \endlatexonly
667	!!
668	!! Calculation of the squared wind shear (required for calculation of the Richardson
669	!! Number)
670	!! \latexonly
671	!! \input{diffucoaero6.tex}
672	!! \endlatexonly
673	!!
674	!! Calculation of the Richardson Number. The Richardson Number is defined as the ratio
675	!! of potential to kinetic energy, or, in the context of atmospheric science, of the
676	!! generation of energy by wind shear against consumption
677	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
678	!! becomes turbulent and vise versa). It is approximated using the expression below:
679	!! \latexonly
680	!! \input{diffucoaero7.tex}
681	!! \endlatexonly
682	!!
683	!! The Richardson Number hence calculated is subject to a minimum value:
684	!! \latexonly
685	!! \input{diffucoaero8.tex}
686	!! \endlatexonly
687	!!
688	!! Computing the drag coefficient. We add the add the height of the vegetation to the
689	!! level height to take into account that the level 'seen' by the vegetation is actually
690	!! the top of the vegetation. Then we we can subtract the displacement height.
691	!! \latexonly
692	!! \input{diffucoaero9.tex}
693	!! \endlatexonly
694	!!
695	!! For the stable case (i.e $R_i$ $\geq$ 0)
696	!! \latexonly
697	!! \input{diffucoaero10.tex}
698	!! \endlatexonly
699	!!
700	!! \latexonly
701	!! \input{diffucoaero11.tex}
702	!! \endlatexonly
703	!!
704	!! For the unstable case (i.e. $R_i$ < 0)
705	!! \latexonly
706	!! \input{diffucoaero12.tex}
707	!! \endlatexonly
708	!!
709	!! \latexonly
710	!! \input{diffucoaero13.tex}
711	!! \endlatexonly
712	!!
713	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
714	!! To prevent this, a minimum limit to the drag coefficient is defined as:
715	!!
716	!! \latexonly
717	!! \input{diffucoaero14.tex}
718	!! \endlatexonly
719	!!
720	!! RECENT CHANGE(S): None
721	!!
722	!! MAIN OUTPUT VARIABLE(S): q_cdrag
723	!!
724	!! REFERENCE(S) :
725	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
726	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
727	!! Circulation Model. Journal of Climate, 6, pp.248-273
728	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
729	!! sur le cycle de l'eau, PhD Thesis, available from:
730	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
731	!! - Jacobson M.Z., Fundamentals of Atmospheric Modeling (2nd Edition), published Cambridge
732	!! University Press, ISBN 0-521-54865-9
733	!!
734	!! FLOWCHART :
735	!! \latexonly
736	!! \includegraphics[scale=0.5]{diffuco_aero_flowchart.png}
737	!! \endlatexonly
738	!! \n
739	!_ ================================================================================================================================
740
741	SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0, roughheight, veget_max, temp_sol, temp_air, &
742	& qsurf, qair, q_cdrag)
743
744	!! 0. Variable and parameter declaration
745
746	!! 0.1 Input variables
747
748	INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size
749	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
750	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
751	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first atmospheric layer (m)
752	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0 !! Surface roughness Length (m)
753	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height (m)
754	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type (-)
755	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Ground temperature (K)
756	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature (K)
757	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific air humidity (kg kg^{-1})
758	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
759
760	!! 0.2 Output variables
761
762	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Product of Surface drag coefficient and wind speed
763	!! (m s^{-1})
764
765	!! 0.3 Modified variables
766
767	!! 0.4 Local variables
768
769	INTEGER(i_std) :: ji, jv
770	REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2
771	REAL(r_std) :: zri, cd_neut, zscf, cd_tmp
772	!_ ================================================================================================================================
773
774	!! 1. Initialisation
775
776	! test if we have to work with q_cdrag or to calcul it
777	DO ji=1,kjpindex
778
779	!! 1a).1 Designation of wind speed
780	!! \latexonly
781	!! \input{diffucoaero1.tex}
782	!! \endlatexonly
783	speed = wind(ji)
784
785	!! 1a).2 Calculation of geopotentiel
786	!! This is the definition of Geopotential height (e.g. Jacobson eqn.4.47, 2005). (required
787	!! for calculation of the Richardson Number)
788	!! \latexonly
789	!! \input{diffucoaero2.tex}
790	!! \endlatexonly
791	zg = zlev(ji) * cte_grav
792
793	!! \latexonly
794	!! \input{diffucoaero3.tex}
795	!! \endlatexonly
796	zdphi = zg/cp_air
797
798	!! 1a).3 Calculation of the virtual air temperature at the surface
799	!! required for calculation of the Richardson Number
800	!! \latexonly
801	!! \input{diffucoaero4.tex}
802	!! \endlatexonly
803	ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji))
804
805	!! 1a).4 Calculation of the virtual surface temperature
806	!! required for calculation of the Richardson Number
807	!! \latexonly
808	!! \input{diffucoaero5.tex}
809	!! \endlatexonly
810	ztvs = temp_sol(ji) * (un + retv * qsurf(ji))
811
812	!! 1a).5 Calculation of the squared wind shear
813	!! required for calculation of the Richardson Number
814	!! \latexonly
815	!! \input{diffucoaero6.tex}
816	!! \endlatexonly
817	zdu2 = MAX(cepdu2,speed**2)
818
819	!! 1a).6 Calculation of the Richardson Number
820	!! The Richardson Number is defined as the ratio of potential to kinetic energy, or, in the
821	!! context of atmospheric science, of the generation of energy by wind shear against consumption
822	!! by static stability and is an indicator of flow stability (i.e. for when laminar flow
823	!! becomes turbulent and vise versa).\n
824	!! It is approximated using the expression below:
825	!! \latexonly
826	!! \input{diffucoaero7.tex}
827	!! \endlatexonly
828	zri = zg * (ztvd - ztvs) / (zdu2 * ztvd)
829
830	!! The Richardson Number hence calculated is subject to a minimum value:
831	!! \latexonly
832	!! \input{diffucoaero8.tex}
833	!! \endlatexonly
834	zri = MAX(MIN(zri,5.),-5.)
835
836	!! 1a).7 Computing the drag coefficient
837	!! We add the add the height of the vegetation to the level height to take into account
838	!! that the level 'seen' by the vegetation is actually the top of the vegetation. Then we
839	!! we can subtract the displacement height.
840	!! \latexonly
841	!! \input{diffucoaero9.tex}
842	!! \endlatexonly
843	cd_neut = (ct_karman / LOG( (zlev(ji) + roughheight(ji)) / z0(ji) )) ** 2
844
845	!! 1a).7.1 - for the stable case (i.e $R_i$ $\geq$ 0)
846	IF (zri .GE. zero) THEN
847
848	!! \latexonly
849	!! \input{diffucoaero10.tex}
850	!! \endlatexonly
851	zscf = SQRT(un + cd * ABS(zri))
852
853	!! \latexonly
854	!! \input{diffucoaero11.tex}
855	!! \endlatexonly
856	cd_tmp=cd_neut/(un + trois * cb * zri * zscf)
857	ELSE
858
859	!! 1a).7.2 - for the unstable case (i.e. $R_i$ < 0)
860	!! \latexonly
861	!! \input{diffucoaero12.tex}
862	!! \endlatexonly
863	zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * &
864	& ((zlev(ji) + roughheight(ji)) / z0(ji))))
865
866	!! \latexonly
867	!! \input{diffucoaero13.tex}
868	!! \endlatexonly
869	cd_tmp=cd_neut * (un - trois * cb * zri * zscf)
870	ENDIF
871
872	!! If the Drag Coefficient becomes too small than the surface may uncouple from the atmosphere.
873	!! To prevent this, a minimum limit to the drag coefficient is defined as:
874
875	!! \latexonly
876	!! \input{diffucoaero14.tex}
877	!! \endlatexonly
878	!!
879	q_cdrag(ji) = MAX(cd_tmp, 1.e-4/MAX(speed,min_wind))
880
881	! In some situations it might be useful to give an upper limit on the cdrag as well.
882	! The line here should then be uncommented.
883	!q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind))
884
885	END DO
886
887	IF (long_print) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done '
888
889	END SUBROUTINE diffuco_aero
890
891
892	!! ================================================================================================================================
893	!! SUBROUTINE : diffuco_snow
894	!!
895	!>\BRIEF This subroutine computes the beta coefficient for snow sublimation.
896	!!
897	!! DESCRIPTION : This routine computes beta coefficient for snow sublimation, which
898	!! integrates the snow on both vegetation and other surface types (e.g. ice, lakes,
899	!! cities etc.) \n
900	!!
901	!! A critical depth of snow (snowcri) is defined to calculate the fraction of each grid-cell
902	!! that is covered with snow (snow/snowcri) while the remaining part is snow-free.
903	!! We also carry out a first calculation of sublimation (subtest) to lower down the beta
904	!! coefficient if necessary (if subtest > snow). This is a predictor-corrector test.
905	!!
906	!! RECENT CHANGE(S): None
907	!!
908	!! MAIN OUTPUT VARIABLE(S): ::vbeta1
909	!!
910	!! REFERENCE(S) :
911	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
912	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
913	!! Circulation Model. Journal of Climate, 6, pp. 248-273
914	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
915	!! sur le cycle de l'eau, PhD Thesis, available from:
916	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
917	!!
918	!! FLOWCHART : None
919	!! \n
920	!_ ================================================================================================================================
921
922	SUBROUTINE diffuco_snow (kjpindex, dtradia, qair, qsatt, rau, u, v,q_cdrag, &
923	& snow, frac_nobio, totfrac_nobio, snow_nobio, vbeta1)
924
925	!! 0. Variable and parameter declaration
926
927	!! 0.1 Input variables
928
929	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
930	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds (s)
931	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
932	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
933	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air density (kg m^{-3})
934	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
935	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
936	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Product of surface drag coefficient and wind speed (s m^{-1})
937	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg m^{-2})
938	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. (-)
939	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice, lakes, cities etc. (-)
940	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice, lakes, cities etc. (-)
941
942	!! 0.2 Output variables
943
944	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation (dimensionless ratio)
945
946	!! 0.3 Modified variables
947
948	!! 0.4 Local variables
949
950	REAL(r_std) :: subtest !! Sublimation for test (kg m^{-2})
951	REAL(r_std) :: zrapp !! Modified factor (ratio)
952	REAL(r_std) :: speed !! Wind speed (m s^{-1})
953	REAL(r_std) :: vbeta1_add !! Beta for sublimation (ratio)
954	INTEGER(i_std) :: ji, jv !! Indices (-)
955	!_ ================================================================================================================================
956
957	!! 1. Calculate beta coefficient for snow sublimation on the vegetation\n
958
959	DO ji=1,kjpindex ! Loop over # pixels - domain size
960
961	! Fraction of mesh that can sublimate snow
962	vbeta1(ji) = (un - totfrac_nobio(ji)) * MAX( MIN(snow(ji)/snowcri,un),zero)
963
964	! Limitation of sublimation in case of snow amounts smaller than the atmospheric demand.
965	speed = MAX(min_wind, wind(ji))
966
967	subtest = dtradia * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * &
968	& ( qsatt(ji) - qair(ji) )
969
970	IF ( subtest .GT. zero ) THEN
971	zrapp = snow(ji) / subtest
972	IF ( zrapp .LT. un ) THEN
973	vbeta1(ji) = vbeta1(ji) * zrapp
974	ENDIF
975	ENDIF
976
977	END DO ! Loop over # pixels - domain size
978
979	!! 2. Add the beta coefficients calculated from other surfaces types (snow on ice,lakes, cities...)
980
981	DO jv = 1, nnobio ! Loop over # other surface types
982	!!$ !
983	!!$ IF ( jv .EQ. iice ) THEN
984	!!$ !
985	!!$ ! Land ice is of course a particular case
986	!!$ !
987	!!$ DO ji=1,kjpindex
988	!!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv)
989	!!$ ENDDO
990	!!$ !
991	!!$ ELSE
992	!
993	DO ji=1,kjpindex ! Loop over # pixels - domain size
994
995	vbeta1_add = frac_nobio(ji,jv) * MAX(MIN(snow_nobio(ji,jv)/snowcri,un), zero)
996
997	! Limitation of sublimation in case of snow amounts smaller than
998	! the atmospheric demand.
999	speed = MAX(min_wind, wind(ji))
1000
1001	!! Limitation of sublimation by the snow accumulated on the ground
1002	!! A first approximation is obtained with the old values of
1003	!! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc)
1004	subtest = dtradia * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * &
1005	& ( qsatt(ji) - qair(ji) )
1006
1007	IF ( subtest .GT. zero ) THEN
1008	zrapp = snow_nobio(ji,jv) / subtest
1009	IF ( zrapp .LT. un ) THEN
1010	vbeta1_add = vbeta1_add * zrapp
1011	ENDIF
1012	ENDIF
1013
1014	vbeta1(ji) = vbeta1(ji) + vbeta1_add
1015
1016	ENDDO ! Loop over # pixels - domain size
1017
1018	!!$ ENDIF
1019
1020	ENDDO ! Loop over # other surface types
1021
1022	IF (long_print) WRITE (numout,*) ' diffuco_snow done '
1023
1024	END SUBROUTINE diffuco_snow
1025
1026
1027	!! ================================================================================================================================
1028	!! SUBROUTINE : diffuco_inter
1029	!!
1030	!>\BRIEF This routine computes the partial beta coefficient
1031	!! for the interception for each type of vegetation
1032	!!
1033	!! DESCRIPTION : We first calculate the dry and wet parts of each PFT (wet part = qsintveg/qsintmax).
1034	!! The former is submitted to transpiration only (vbeta3 coefficient, calculated in
1035	!! diffuco_trans or diffuco_trans_co2), while the latter is first submitted to interception loss
1036	!! (vbeta2 coefficient) and then to transpiration once all the intercepted water has been evaporated
1037	!! (vbeta23 coefficient). Interception loss is also submitted to a predictor-corrector test,
1038	!! as for snow sublimation. \n
1039	!!
1040	!! \latexonly
1041	!! \input{diffucointer1.tex}
1042	!! \endlatexonly
1043	!! Calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1044	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1045	!!
1046	!! \latexonly
1047	!! \input{diffucointer2.tex}
1048	!! \endlatexonly
1049	!!
1050	!! Calculation of $\beta_3$, the canopy transpiration resistance
1051	!! \latexonly
1052	!! \input{diffucointer3.tex}
1053	!! \endlatexonly
1054	!!
1055	!! We here determine the limitation of interception loss by the water stored on the leaf.
1056	!! A first approximation of interception loss is obtained using the old values of
1057	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1058	!! \latexonly
1059	!! \input{diffucointer4.tex}
1060	!! \endlatexonly
1061	!!
1062	!! \latexonly
1063	!! \input{diffucointer5.tex}
1064	!! \endlatexonly
1065	!!
1066	!! \latexonly
1067	!! \input{diffucointer6.tex}
1068	!! \endlatexonly
1069	!!
1070	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1071	!! 'zqsvegrap'.
1072	!! \latexonly
1073	!! \input{diffucointer7.tex}
1074	!! \endlatexonly
1075	!!
1076	!! RECENT CHANGE(S): None
1077	!!
1078	!! MAIN OUTPUT VARIABLE(S): ::vbeta2, ::vbeta23
1079	!!
1080	!! REFERENCE(S) :
1081	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1082	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1083	!! Circulation Model. Journal of Climate, 6, pp. 248-273
1084	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1085	!! sur le cycle de l'eau, PhD Thesis, available from:
1086	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1087	!! - Perrier, A, 1975. Etude physique de l'Ã©vaporation dans les conditions naturelles. Annales
1088	!! Agronomiques, 26(1-18): pp. 105-123, pp. 229-243
1089	!!
1090	!! FLOWCHART : None
1091	!! \n
1092	!_ ================================================================================================================================
1093
1094	SUBROUTINE diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
1095	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
1096
1097	! Nathalie - Juin 2006 - Introduction de vbeta23
1098	! SUBROUTINE diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
1099	! & qsintveg, qsintmax, rstruct, vbeta2)
1100
1101	!! 0 Variable and parameter declaration
1102
1103	!! 0.1 Input variables
1104
1105	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1106	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
1107	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
1108	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity (kg kg^{-1})
1109	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
1110	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
1111	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
1112	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Product of Surface drag coefficient and wind
1113	!! speed (m s^{-1})
1114	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type (fraction)
1115	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
1116	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
1117	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! architectural resistance (s m^{-1})
1118
1119	!! 0.2 Output variables
1120
1121	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss (-)
1122	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23 !! Beta for fraction of wetted foliage that will
1123	!! transpire (-)
1124
1125	!! 0.4 Local variables
1126
1127	INTEGER(i_std) :: ji, jv !! (-), (-)
1128	REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed !!
1129	!_ ================================================================================================================================
1130
1131	!! 1. Initialize
1132
1133	vbeta2(:,:) = zero
1134	vbeta23(:,:) = zero
1135
1136	!! 2. The beta coefficient for interception by vegetation.
1137
1138	DO jv = 1,nvm
1139
1140	DO ji=1,kjpindex
1141
1142	IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN
1143
1144	zqsvegrap = zero
1145	IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN
1146
1147	!! \latexonly
1148	!! \input{diffucointer1.tex}
1149	!! \endlatexonly
1150	!!
1151	!! We calculate the wet fraction of vegetation as the ration between the intercepted water and the maximum water
1152	!! on the vegetation. This ratio defines the wet portion of vegetation that will be submitted to interception loss.
1153	!!
1154	zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
1155	END IF
1156	!! \latexonly
1157	!! \input{diffucointer2.tex}
1158	!! \endlatexonly
1159	speed = MAX(min_wind, wind(ji))
1160
1161	!! Calculation of $\beta_3$, the canopy transpiration resistance
1162	!! \latexonly
1163	!! \input{diffucointer3.tex}
1164	!! \endlatexonly
1165	vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv)))
1166
1167	!! We here determine the limitation of interception loss by the water stored on the leaf.
1168	!! A first approximation of interception loss is obtained using the old values of
1169	!! qair and qsol_sat, which are functions of temp-sol and pb. (see call of 'qsatcalc')
1170	!! \latexonly
1171	!! \input{diffucointer4.tex}
1172	!! \endlatexonly
1173	ziltest = dtradia * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * &
1174	& ( qsatt(ji) - qair(ji) )
1175
1176	IF ( ziltest .GT. zero ) THEN
1177
1178	!! \latexonly
1179	!! \input{diffucointer5.tex}
1180	!! \endlatexonly
1181	zrapp = qsintveg(ji,jv) / ziltest
1182	IF ( zrapp .LT. un ) THEN
1183
1184	!! \latexonly
1185	!! \input{diffucointer6.tex}
1186	!! \endlatexonly
1187	!!
1188	!! Once the whole water stored on foliage has evaporated, transpiration can take place on the fraction
1189	!! 'zqsvegrap'.
1190	vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, zero)
1191
1192	!! \latexonly
1193	!! \input{diffucointer7.tex}
1194	!! \endlatexonly
1195	vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp
1196	ENDIF
1197	ENDIF
1198	END IF
1199
1200	END DO
1201
1202	END DO
1203
1204	IF (long_print) WRITE (numout,*) ' diffuco_inter done '
1205
1206	END SUBROUTINE diffuco_inter
1207
1208
1209	!! ==============================================================================================================================
1210	!! SUBROUTINE : diffuco_bare
1211	!!
1212	!>\BRIEF This routine computes the partial beta coefficient corresponding to
1213	!! bare soil
1214	!!
1215	!! DESCRIPTION : Bare soil evaporation is either limited by a soil resistance
1216	!! (rsol) that is calculated in hydrolc.f90, when Choisnel hydrology is used or
1217	!! submitted to a maximum possible flow (evap_bare_lim) if the 11-layer hydrology is used.\n
1218	!!
1219	!! Calculation of wind speed
1220	!! \latexonly
1221	!! \input{diffucobare1.tex}
1222	!! \endlatexonly
1223	!!
1224	!! The calculation of $\beta_4$
1225	!! \latexonly
1226	!! \input{diffucobare2.tex}
1227	!! \endlatexonly
1228	!!
1229	!! RECENT CHANGE(S): None
1230	!!
1231	!! MAIN OUTPUT VARIABLE(S): ::vbeta4
1232	!!
1233	!! REFERENCE(S) :
1234	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1235	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1236	!! Circulation Model. Journal of Climate, 6, pp.248-273
1237	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1238	!! sur le cycle de l'eau, PhD Thesis, available from:
1239	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1240	!!
1241	!! FLOWCHART : None
1242	!! \n
1243	!_ ================================================================================================================================
1244
1245	SUBROUTINE diffuco_bare (kjpindex, dtradia, u, v, q_cdrag, rsol, evap_bare_lim, evapot, humrel, veget, vbeta4)
1246
1247	!! 0. Variable and parameter declaration
1248
1249	!! 0.1 Input variables
1250	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1251	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
1252	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
1253	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
1254	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Product of Surface drag coefficient and wind speed
1255	!! (m s^{-1})
1256	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rsol !! resistance for bare soil evaporation (s m^{-1})
1257	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! limiting factor for bare soil evaporation when the
1258	!! 11-layer hydrology is used (-)
1259	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation (??mm day^{-1}??)
1260	!! NdN Variable that does not seem to be used at
1261	!! all in diffuco!!
1262	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
1263	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction (-)
1264
1265	!! 0.2 Output variables
1266
1267	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation (-)
1268
1269	!! 0.3 Modified variables
1270
1271	!! 0.4 Local variables
1272
1273	INTEGER(i_std) :: ji, jv
1274	REAL(r_std) :: speed !! Surface wind speed (m s^{-1})
1275	!_ ================================================================================================================================
1276
1277	!! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil
1278
1279	IF ( .NOT. control%hydrol_cwrr ) THEN
1280	DO ji = 1, kjpindex
1281
1282	vbeta4(ji) = zero
1283
1284	! Calculation of the soil resistance and the beta ($\beta_4$) for bare soil.
1285	! note: veget ( ,1) contains the fraction of bare soil
1286	! see hydrol module
1287	IF (veget(ji,1) .GE. min_sechiba) THEN
1288
1289	!! \latexonly
1290	!! \input{diffucobare1.tex}
1291	!! \endlatexonly
1292	speed = MAX(min_wind, wind(ji))
1293
1294	! Correction Nathalie de Noblet - le 27 Mars 2006
1295	! Selon recommandation de Frederic Hourdin: supprimer humrel dans formulation vbeta4
1296	! vbeta4(ji) = veget(ji,1) humrel(ji,1) (un / (un + speed * q_cdrag(ji) * rsol(ji)))
1297	! Nathalie - le 28 mars 2006 - vbeta4 n'etait pas calcule en fonction de
1298	! rsol mais de rsol_cste * hdry! Dans ce cas inutile de calculer rsol(ji)!!
1299
1300	!! The calculation of $\beta_4$
1301	!! \latexonly
1302	!! \input{diffucobare2.tex}
1303	!! \endlatexonly
1304	vbeta4(ji) = veget(ji,1) * (un / (un + speed * q_cdrag(ji) * rsol(ji)))
1305
1306	ENDIF
1307
1308	END DO
1309	ELSE
1310	DO ji = 1, kjpindex
1311
1312	!! \latexonly
1313	!! \input{diffucobare3.tex}
1314	!! \endlatexonly
1315	vbeta4(ji) = evap_bare_lim(ji)
1316	END DO
1317	ENDIF
1318
1319	IF (long_print) WRITE (numout,*) ' diffuco_bare done '
1320
1321	END SUBROUTINE diffuco_bare
1322
1323
1324	!! ================================================================================================================================
1325	!! SUBROUTINE : diffuco_trans
1326	!!
1327	!>\BRIEF This routine computes the partial beta coefficient
1328	!! corresponding to transpiration for each vegetation type.
1329	!!
1330	!! DESCRIPTION : Beta coefficients for transpiration are calculated
1331	!! here using Jarvis formulation for stomatal resistance and
1332	!! the structural resistance to represent the vertical gradient of
1333	!! transpiration within the canopy. \n
1334	!!
1335	!! The Jarvis formulation as used here is derived by Lohanner et al. (1980) from Jarvis (1976). This formulation is
1336	!! semi-empirical: \n
1337	!!
1338	!! \latexonly
1339	!! \input{diffucotrans4.tex}
1340	!! \endlatexonly
1341	!! \n
1342	!!
1343	!! where in this expression LAI is the single sided Leaf Area Index, R_{new}^{SW} the net shortwave radiation,
1344	!! R_{SO} the half light saturation factor, \delta c the water vapour concentration deficit, and a, k_0 and \lambda
1345	!! are all parameters that are derived from extensive measurement of surface layer vegetation. \n
1346	!!
1347	!! Structural resistance (or architectural resistance) is a function of vegetation type and is assigned based on the
1348	!! particular Plant Functional Type (PFT) in question. The range of values for the structural resistance are listed
1349	!! in the module 'pft_parameters', and are described in de Noblet-DucoudrÃ© et al (1993). \n
1350	!!
1351	!! vbetaco2 is here set to zero as this way to compute canopy resistances is only used
1352	!! without STOMATE, and there is therefore no photosynthesis. \n
1353	!!
1354	!! Moisture concentration at the leaf level.
1355	!! \latexonly
1356	!! \input{diffucotrans1.tex}
1357	!! \endlatexonly
1358	!!
1359	!! Calulation of the beta coefficient for vegetation transpiration, beta_3.
1360	!! \latexonly
1361	!! \input{diffucotrans2.tex}
1362	!! \endlatexonly
1363	!! \latexonly
1364	!! \input{diffucotrans3.tex}
1365	!! \endlatexonly
1366	!!
1367	!! \latexonly
1368	!! \input{diffucotrans4.tex}
1369	!! \endlatexonly
1370	!!
1371	!! This is the formulation for beta_3.
1372	!! \latexonly
1373	!! \input{diffucotrans5.tex}
1374	!! \endlatexonly
1375	!!
1376	!! \latexonly
1377	!! \input{diffucotrans6.tex}
1378	!! \endlatexonly
1379	!!
1380	!! RECENT CHANGE(S): None
1381	!!
1382	!! MAIN OUTPUT VARIABLE(S): ::vbeta3, ::rveget, ::cimean and ::vbetaco2
1383	!!
1384	!! REFERENCE(S) :
1385	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
1386	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
1387	!! Circulation Model. Journal of Climate, 6, pp.248-273
1388	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
1389	!! sur le cycle de l'eau, PhD Thesis, available from:
1390	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
1391	!! - Jarvis, PG, 1976. The interpretation of the variations in leaf water potential and stomatal
1392	!! conductance found in canopies in the fields. Philosophical Transactions of the Royal Society of
1393	!! London, Series B, 273, pp. 593-610
1394	!! - Lohammer T, Larsson S, Linder S & Falk SO, 1980. Simulation models of gaseous exchange in Scotch
1395	!! pine. Structure and function of Northern Coniferous Forest, Ecological Bulletin, 32, pp. 505-523
1396	!!
1397	!! FLOWCHART : None
1398	!! \n
1399	!_ ================================================================================================================================
1400
1401	SUBROUTINE diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
1402	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
1403
1404	! Nathalie - Juin 2006 - introduction de vbeta23
1405	! SUBROUTINE diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
1406	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
1407
1408
1409	!! 0. Variable and parameter declaration
1410
1411	!! 0.1 Input variables
1412
1413	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
1414	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
1415	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Short wave net flux at surface (W m^{-2})
1416	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
1417	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (Pa)
1418	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
1419	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
1420	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
1421	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind speed (m s^{-1})
1422	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Product of Surface drag coefficient and wind speed
1423	!! (m s^{-1})
1424	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (within range 0 to 1)
1425	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation (-)
1426	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. vegetation fraction (LAI->infty) (fraction)
1427	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
1428	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception (kg m^{-2})
1429	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
1430	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! Structural resistance (s m^{-1})
1431	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for wetted foliage fraction that will transpire (-)
1432
1433	!! 0.2 Output variables
1434
1435	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration (-)
1436	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Stomatal resistance of the whole canopy (s m^{-1})
1437	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE: mean intercellular ci (see enerbil)
1438	!! (\mumole m^{-2} s^{-1})
1439	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE: Beta for CO2 (-)
1440
1441	!! 0.3 Modified variables
1442
1443	!! 0.4 Local variables
1444
1445	INTEGER(i_std) :: ji, jv
1446	REAL(r_std) :: speed
1447	REAL(r_std), DIMENSION(kjpindex) :: zdefconc, zqsvegrap
1448	REAL(r_std), DIMENSION(kjpindex) :: qsatt
1449	!_ ================================================================================================================================
1450
1451
1452	!! 1. Moisture concentration at the leaf level.
1453
1454	CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
1455
1456	!! \latexonly
1457	!! \input{diffucotrans1.tex}
1458	!! \endlatexonly
1459	zdefconc(:) = rau(:) * MAX( qsatt(:) - qair(:), zero )
1460
1461
1462	!! 2. Calulation of the beta coefficient for vegetation transpiration, beta_3.
1463
1464	DO jv = 1,nvm
1465
1466	rveget(:,jv) = undef_sechiba
1467	vbeta3(:,jv) = zero
1468
1469	zqsvegrap(:) = zero
1470
1471	DO ji = 1, kjpindex
1472
1473	!! \latexonly
1474	!! \input{diffucotrans2.tex}
1475	!! \endlatexonly
1476	speed = MAX(min_wind, wind(ji))
1477
1478	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
1479
1480	!! \latexonly
1481	!! \input{diffucotrans3.tex}
1482	!! \endlatexonly
1483	zqsvegrap(ji) = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
1484	ENDIF
1485
1486	IF ( ( veget(ji,jv)*lai(ji,jv) .GT. min_sechiba ) .AND. &
1487	( kzero(jv) .GT. min_sechiba ) .AND. &
1488	( swnet(ji) .GT. min_sechiba ) ) THEN
1489
1490	!! \latexonly
1491	!! \input{diffucotrans4.tex}
1492	!! \endlatexonly
1493	rveget(ji,jv) = (( swnet(ji) + rayt_cste ) / swnet(ji) ) &
1494	* ((defc_plus + (defc_mult * zdefconc(ji) )) / kzero(jv)) * (un / lai(ji,jv))
1495
1496	! Corrections Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
1497	! Introduction d'un potentiometre (rveg_pft) pour regler la somme rveg+rstruct
1498	! vbeta3(ji,jv) = veget(ji,jv) * (un - zqsvegrap(ji)) * humrel(ji,jv) * &
1499	! (un / (un + speed * q_cdrag(ji) * (rveget(ji,jv) + rstruct(ji,jv))))
1500
1501	!! This is the formulation for $beta_3$.
1502	!! \latexonly
1503	!! \input{diffucotrans5.tex}
1504	!! \endlatexonly
1505	vbeta3(ji,jv) = veget(ji,jv) * (un - zqsvegrap(ji)) * humrel(ji,jv) * &
1506	(un / (un + speed * q_cdrag(ji) * (rveg_pft(jv)*(rveget(ji,jv) + rstruct(ji,jv)))))
1507
1508	! Fin ajout Nathalie
1509	! Ajout Nathalie - Juin 2006
1510
1511	!! \latexonly
1512	!! \input{diffucotrans6.tex}
1513	!! \endlatexonly
1514	vbeta3(ji,jv) = vbeta3(ji,jv) + MIN( vbeta23(ji,jv), &
1515	veget(ji,jv) * zqsvegrap(ji) * humrel(ji,jv) * &
1516	(un / (un + speed * q_cdrag(ji) * (rveg_pft(jv)*(rveget(ji,jv) + rstruct(ji,jv))))))
1517	! Fin ajout Nathalie
1518
1519	ENDIF
1520
1521	ENDDO
1522
1523	ENDDO
1524
1525	! STOMATE
1526	cimean(:,:) = zero
1527	vbetaco2(:,:) = zero
1528
1529	IF (long_print) WRITE (numout,*) ' diffuco_trans done '
1530
1531	END SUBROUTINE diffuco_trans
1532
1533
1534	!! ==============================================================================================================================
1535	!! SUBROUTINE : diffuco_trans_co2
1536	!!
1537	!>\BRIEF This subroutine computes carbon assimilation and stomatal
1538	!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987).
1539	!!
1540	!! DESCRIPTION :\n
1541	!! *** General:\n
1542	!! The equations are different depending on the photosynthesis mode (C3 versus C4).\n
1543	!! Assimilation and conductance are computed over 20 levels of LAI and then
1544	!! integrated at the canopy level.\n
1545	!! This routine also computes partial beta coefficient: transpiration for each
1546	!! type of vegetation.\n
1547	!! There is a main loop on the PFTs, then inner loops on the points where
1548	!! assimilation has to be calculated.\n
1549	!! This subroutine is called by diffuco_main only if photosynthesis is activated
1550	!! for sechiba (flag STOMATE_OK_CO2=TRUE), otherwise diffuco_trans is called.\n
1551	!! This subroutine is called at each sechiba time step by sechiba_main.\n
1552	!! *** Details:
1553	!! - Integration at the canopy level\n
1554	!! \latexonly
1555	!! \input{diffuco_trans_co2_1.1.tex}
1556	!! \endlatexonly
1557	!! - Light''s extinction \n
1558	!! The available light follows a simple Beer extinction law.
1559	!! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.\n
1560	!! \latexonly
1561	!! \input{diffuco_trans_co2_1.2.tex}
1562	!! \endlatexonly
1563	!! - Estimation of relative humidity of air (for calculation of the stomatal conductance)\n
1564	!! \latexonly
1565	!! \input{diffuco_trans_co2_1.3.tex}
1566	!! \endlatexonly
1567	!! - Calculation of the water limitation factor\n
1568	!! \latexonly
1569	!! \input{diffuco_trans_co2_2.1.tex}
1570	!! \endlatexonly
1571	!! - Calculation of temperature dependent parameters for C4 plants\n
1572	!! \latexonly
1573	!! \input{diffuco_trans_co2_2.2.tex}
1574	!! \endlatexonly
1575	!! - Calculation of temperature dependent parameters for C3 plants\n
1576	!! \latexonly
1577	!! \input{diffuco_trans_co2_2.3.tex}
1578	!! \endlatexonly
1579	!! - Vmax scaling\n
1580	!! Vmax is scaled into the canopy due to reduction of nitrogen
1581	!! (Johnson and Thornley,1984).\n
1582	!! \latexonly
1583	!! \input{diffuco_trans_co2_2.4.1.tex}
1584	!! \endlatexonly
1585	!! - Assimilation for C4 plants (Collatz et al., 1992)\n
1586	!! \latexonly
1587	!! \input{diffuco_trans_co2_2.4.2.tex}
1588	!! \endlatexonly
1589	!! - Assimilation for C3 plants (Farqhuar et al., 1980)\n
1590	!! \latexonly
1591	!! \input{diffuco_trans_co2_2.4.3.tex}
1592	!! \endlatexonly
1593	!! - Estimation of the stomatal conductance (Ball et al., 1987)\n
1594	!! \latexonly
1595	!! \input{diffuco_trans_co2_2.4.4.tex}
1596	!! \endlatexonly
1597	!!
1598	!! RECENT CHANGE(S): N. de Noblet 2006/06
1599	!! - addition of q2m and t2m as input parameters for the
1600	!! calculation of Rveget
1601	!! - introduction of vbeta23
1602	!!
1603	!! MAIN OUTPUT VARIABLE(S): beta coefficients, resistances, CO2 intercellular
1604	!! concentration
1605	!!
1606	!! REFERENCE(S) :
1607	!! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal
1608	!! conductance and its contribution to the control of photosynthesis under
1609	!! different environmental conditions, Prog. Photosynthesis, 4, 221â 224.
1610	!! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis
1611	!! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol.,
1612	!! 19, 519â538.
1613	!! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of
1614	!! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78â90.
1615	!! - Johnson, I. R., and J. Thornley (1984), A model of instantaneous and daily
1616	!! canopy photosynthesis, J Theor. Biol., 107, 531 545
1617	!! - McMurtrie, R.E., Rook, D.A. and Kelliher, F.M., 1990. Modelling the yield of Pinus radiata on a
1618	!! site limited by water and nitrogen. For. Ecol. Manage., 30: 381-413
1619	!!
1620	!! FLOWCHART : None
1621	!! \n
1622	!_ ================================================================================================================================
1623
1624	SUBROUTINE diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, &
1625	assim_param, ccanopy, &
1626	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
1627
1628	!
1629	!! 0. Variable and parameter declaration
1630	!
1631
1632	!
1633	!! 0.1 Input variables
1634	!
1635	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1636	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
1637	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Downwelling short wave flux
1638	!! @tex ($W m^{-2}$) @endtex
1639	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature (K)
1640	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (Pa)
1641	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
1642	!! @tex ($kg kg^{-1}$) @endtex
1643	! N. de Noblet - 2006/06 - addition of q2m and t2m
1644	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity
1645	!! @tex ($kg kg^{-1}$) @endtex
1646	! In off-line mode q2m and qair are the same.
1647	! In off-line mode t2m and temp_air are the same.
1648	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature (K)
1649	! N. de Noblet - 2006/06 - addition of q2m and t2m
1650	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! air density @tex ($kg m^{-3}$) @endtex
1651	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
1652	!! @tex ($m s^{-1}$) @endtex
1653	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
1654	!! @tex ($m s^{-1}$) @endtex
1655	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
1656	!! @tex ($m s^{-1}$) @endtex
1657	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param !! min+max+opt temps (K), vcmax, vjmax for
1658	!! photosynthesis
1659	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1660	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy (ppm)
1661	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress (0-1,unitless)
1662	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Coverage fraction of vegetation for each PFT
1663	!! depending on LAI (0-1, unitless)
1664	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum vegetation fraction of each PFT inside
1665	!! the grid box (0-1, unitless)
1666	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index @tex ($m^2 m^{-2}$) @endtex
1667	!! @tex ($m s^{-1}$) @endtex
1668	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
1669	!! @tex ($kg m^{-2}$) @endtex
1670	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
1671	!! @tex ($kg m^{-2}$) @endtex
1672	! N. de Noblet - 2006/06
1673	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23 !! Beta for fraction of wetted foliage that will
1674	!! transpire (unitless)
1675	! N. de Noblet - 2006/06
1676
1677	!
1678	!! 0.2 Output variables
1679	!
1680	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration (unitless)
1681	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! stomatal resistance of vegetation
1682	!! @tex ($s m^{-1}$) @endtex
1683	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! structural resistance @tex ($s m^{-1}$) @endtex
1684	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! mean intercellular CO2 concentration
1685	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1686	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! beta for CO2 transfer (unitless)
1687
1688	!
1689	!! 0.3 Modified variables
1690	!
1691
1692	!
1693	!! 0.4 Local variables
1694	!
1695	REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax !! maximum rate of carboxylation
1696	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1697	REAL(r_std),DIMENSION (kjpindex,nvm) :: vjmax !! maximum rate of Rubisco regeneration
1698	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1699	REAL(r_std),DIMENSION (kjpindex,nvm) :: tmin !! minimum temperature for photosynthesis (K)
1700	REAL(r_std),DIMENSION (kjpindex,nvm) :: topt !! optimal temperature for photosynthesis (K)
1701	REAL(r_std),DIMENSION (kjpindex,nvm) :: tmax !! maximum temperature for photosynthesis (K)
1702	INTEGER(i_std) :: ji, jv, jl !! indices (unitless)
1703	REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest !! intercellular CO2 concentration at the lowest
1704	!! LAI level (ppm)
1705	INTEGER(i_std), DIMENSION(kjpindex) :: ilai !! counter for loops on LAI levels (unitless)
1706	REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap !! relative water quantity in the water
1707	!! interception reservoir (0-1,unitless)
1708	REAL(r_std) :: speed !! wind speed @tex ($m s^{-1}$) @endtex
1709	! Assimilation
1710	LOGICAL, DIMENSION(kjpindex) :: assimilate !! where assimilation is to be calculated
1711	!! (unitless)
1712	LOGICAL, DIMENSION(kjpindex) :: calculate !! where assimilation is to be calculated for
1713	!! in the PFTs loop (unitless)
1714	INTEGER(i_std) :: nic,inic,icinic !! counter/indices (unitless)
1715	INTEGER(i_std), DIMENSION(kjpindex) :: index_calc !! index (unitless)
1716	INTEGER(i_std) :: nia,inia,nina,inina,iainia !! counter/indices (unitless)
1717	INTEGER(i_std), DIMENSION(kjpindex) :: index_assi,index_non_assi !! indices (unitless)
1718	REAL(r_std), DIMENSION(kjpindex) :: vc2 !! rate of carboxylation (at a specific LAI level)
1719	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1720	REAL(r_std), DIMENSION(kjpindex) :: vj2 !! rate of Rubisco regeneration (at a specific LAI
1721	!! level) @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1722	REAL(r_std), DIMENSION(kjpindex) :: assimi !! assimilation (at a specific LAI level)
1723	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1724	REAL(r_std) :: x_1,x_2,x_3,x_4,x_5,x_6 !! terms of second order polynomials
1725	!! (temporary variables)
1726	REAL(r_std), DIMENSION(kjpindex) :: gstop !! stomatal conductance at topmost level
1727	!! @tex ($m s^{-1}$) @endtex
1728	REAL(r_std), DIMENSION(kjpindex) :: gs !! stomatal conductance
1729	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1730	REAL(r_std), DIMENSION(kjpindex) :: Kc !! Michaelis-Menten constant for the rubisco
1731	!! enzyme catalytic activity for CO2 (ppm)
1732	REAL(r_std), DIMENSION(kjpindex) :: Ko !! Michaelis-Menten constant for the rubisco
1733	!! enzyme catalytic activity for O2 (ppm)
1734	REAL(r_std), DIMENSION(kjpindex) :: CP !! CO2 compensation point (ppm)
1735	REAL(r_std), DIMENSION(kjpindex) :: vc !! rate of carboxylation
1736	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1737	REAL(r_std), DIMENSION(kjpindex) :: vj !! rate of Rubisco regeneration
1738	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1739	REAL(r_std), DIMENSION(kjpindex) :: kt !! temperature-dependent pseudo-first-order rate
1740	!! constant of assimilation response to CO2 (C4)
1741	!! (??)
1742	REAL(r_std), DIMENSION(kjpindex) :: rt !! temperature-dependent non-photosynthetic
1743	!! respiration (C4) @tex ($\mu mol $) @endtex ??
1744	REAL(r_std), DIMENSION(kjpindex) :: air_relhum !! air relative humidity at 2m
1745	!! @tex ($kg kg^{-1}$) @endtex
1746	REAL(r_std), DIMENSION(kjpindex) :: water_lim !! water limitation factor (0-1,unitless)
1747	REAL(r_std), DIMENSION(kjpindex) :: temp_lim !! temperature limitation factor (0-1,unitles)
1748	REAL(r_std), DIMENSION(kjpindex) :: gstot !! total stomatal conductance
1749	!! @tex ($m s^{-1}$) @endtex
1750	REAL(r_std), DIMENSION(kjpindex) :: assimtot !! total assimilation
1751	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1752	REAL(r_std), DIMENSION(kjpindex) :: leaf_gs_top !! leaf stomatal conductance at topmost level
1753	!! @tex ($\mu mol m^{-2} s^{-1}$) @endtex
1754	REAL(r_std), DIMENSION(nlai+1) :: laitab !! tabulated LAI steps @tex ($m^2 m^{-2}$) @endtex
1755	REAL(r_std), DIMENSION(kjpindex) :: qsatt !! surface saturated humidity at 2m (??)
1756	!! @tex ($g g^{-1}$) @endtex
1757	REAL(r_std), DIMENSION(nvm,nlai) :: light !! fraction of light that gets through upper LAI
1758	!! levels (0-1,unitless)
1759	REAL(r_std), DIMENSION(kjpindex) :: ci_gs !! intermediate calculus variable
1760	!! for vectorization
1761	REAL(r_std) :: cresist !! coefficient for resistances (??)
1762
1763	! PARAMETER VALUES
1764	REAL(r_std), PARAMETER :: laimax = 12. !! maximum LAI @tex ($m^2 m^{-2}$) @endtex
1765	REAL(r_std), PARAMETER :: xc4_1 = .83
1766	REAL(r_std), PARAMETER :: xc4_2 = .93
1767
1768	! @defgroup Photosynthesis Photosynthesis
1769	! @{
1770	! 1. Preliminary calculations\n
1771	!
1772	! 1.1 Calculate LAI steps\n
1773	! The integration at the canopy level is done over nlai fixed levels.
1774	!! \latexonly
1775	!! \input{diffuco_trans_co2_1.1.tex}
1776	!! \endlatexonly
1777	! @}
1778	! @codeinc
1779	DO jl = 1, nlai+1
1780	laitab(jl) = laimax(EXP(.15REAL(jl-1,r_std))-1.)/(EXP(.15*REAL(nlai,r_std))-un)
1781	ENDDO
1782	! @endcodeinc
1783
1784	! @addtogroup Photosynthesis
1785	! @{
1786	!
1787	! 1.2 Calculate light fraction for each LAI step\n
1788	! The available light follows a simple Beer extinction law.
1789	! The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.
1790	!! \latexonly
1791	!! \input{diffuco_trans_co2_1.2.tex}
1792	!! \endlatexonly
1793	! @}
1794	! @codeinc
1795	DO jl = 1, nlai
1796	DO jv = 1, nvm
1797	light(jv,jl) = exp( -ext_coef(jv)*laitab(jl) )
1798	ENDDO
1799	ENDDO
1800	! @endcodeinc
1801	!
1802	! Photosynthesis parameters
1803	!
1804	! temperatures in K
1805	!
1806	tmin(:,:) = assim_param(:,:,itmin)
1807	tmax(:,:) = assim_param(:,:,itmax)
1808	topt(:,:) = assim_param(:,:,itopt)
1809	!
1810	vcmax(:,:) = assim_param(:,:,ivcmax)
1811	vjmax(:,:) = assim_param(:,:,ivjmax)
1812
1813	! @addtogroup Photosynthesis
1814	! @{
1815	!
1816	! 1.3 Estimate relative humidity of air (for calculation of the stomatal conductance).\n
1817	!! \latexonly
1818	!! \input{diffuco_trans_co2_1.3.tex}
1819	!! \endlatexonly
1820	! @}
1821	!
1822	! N. de Noblet - 2006/06 - We use q2m/t2m instead of qair.
1823	! CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
1824	! air_relhum(:) = &
1825	! ( qair(:) * pb(:) / (0.622+qair(:)*0.378) ) / &
1826	! ( qsatt(:)pb(:) / (0.622+qsatt(:)0.378 ) )
1827	! @codeinc
1828	CALL qsatcalc (kjpindex, t2m, pb, qsatt)
1829	air_relhum(:) = &
1830	( q2m(:) * pb(:) / (0.622+q2m(:)*0.378) ) / &
1831	( qsatt(:)pb(:) / (0.622+qsatt(:)0.378 ) )
1832	! @endcodeinc
1833	! N. de Noblet - 2006/06
1834	!
1835	! @addtogroup Photosynthesis
1836	! @{
1837	! 2. Loop over vegetation types\n
1838	! @}
1839	!
1840	DO jv = 1,nvm
1841	!
1842	! @addtogroup Photosynthesis
1843	! @{
1844	!
1845	! 2.1 Initializations\n
1846	!! \latexonly
1847	!! \input{diffuco_trans_co2_2.1.tex}
1848	!! \endlatexonly
1849	! @}
1850	!
1851	! beta coefficient for vegetation transpiration
1852	!
1853	rstruct(:,jv) = rstruct_const(jv)
1854	rveget(:,jv) = undef_sechiba
1855	!
1856	vbeta3(:,jv) = zero
1857	vbetaco2(:,jv) = zero
1858	!
1859	cimean(:,jv) = ccanopy(:)
1860	!
1861	!! mask that contains points where there is photosynthesis
1862	!! For the sake of vectorisation [DISPENSABLE], computations are done only for convenient points.
1863	!! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not
1864	!! calculated (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation,
1865	!! temperature and relative humidity).
1866	!! For the points where assimilation is not calculated, variables are initialized to specific values.
1867	!! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
1868	!! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes
1869	!! the nina points with no assimilation.
1870	nia=0
1871	nina=0
1872	!
1873	DO ji=1,kjpindex
1874	!
1875	IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. &
1876	( veget_max(ji,jv) .GT. 1.E-8 ) ) THEN
1877	IF ( ( veget(ji,jv) .GT. 1.E-8 ) .AND. &
1878	( swdown(ji) .GT. min_sechiba ) .AND. &
1879	( temp_air(ji) .GT. tmin(ji,jv) ) .AND. &
1880	( temp_air(ji) .LT. tmax(ji,jv) ) .AND. &
1881	( humrel(ji,jv) .GT. min_sechiba ) ) THEN
1882	!
1883	assimilate(ji) = .TRUE.
1884	nia=nia+1
1885	index_assi(nia)=ji
1886	!
1887	ELSE
1888	!
1889	assimilate(ji) = .FALSE.
1890	nina=nina+1
1891	index_non_assi(nina)=ji
1892	!
1893	ENDIF
1894	ELSE
1895	!
1896	assimilate(ji) = .FALSE.
1897	nina=nina+1
1898	index_non_assi(nina)=ji
1899	!
1900	ENDIF
1901	!
1902	ENDDO
1903	!
1904	gstot(:) = zero
1905	assimtot(:) = zero
1906	!
1907	zqsvegrap(:) = zero
1908	WHERE (qsintmax(:,jv) .GT. min_sechiba)
1909	!! relative water quantity in the water interception reservoir
1910	zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
1911	ENDWHERE
1912	!
1913	!! Calculates the water limitation factor.
1914	WHERE ( assimilate(:) )
1915	water_lim(:) = MIN( 2.*humrel(:,jv), un )
1916	ENDWHERE
1917	! give a default value of ci for all pixel that do not assimilate
1918	DO jl=1,nlai
1919	DO inina=1,nina
1920	leaf_ci(index_non_assi(inina),jv,jl) = ccanopy(index_non_assi(inina)) * .667_r_std
1921	ENDDO
1922	ENDDO
1923	!
1924	ilai(:) = 1
1925	!
1926	! Here is calculated photosynthesis (Farqhuar et al., 1980)
1927	! and stomatal conductance (Ball et al., 1987)
1928	!
1929	! Calculates temperature dependent parameters.
1930	!
1931	IF ( is_c4(jv) ) THEN
1932	!
1933	! @addtogroup Photosynthesis
1934	! @{
1935	!
1936	! 2.2 Calculates temperature dependent parameters for C4 plants.\n
1937	!! \latexonly
1938	!! \input{diffuco_trans_co2_2.2.tex}
1939	!! \endlatexonly
1940	! @}
1941	!
1942	IF (nia .GT. 0) then
1943	DO inia=1,nia
1944	!
1945	! @codeinc
1946	x_1 = 0.177 * EXP( 0.069*(temp_air(index_assi(inia))-tp_00) )
1947	! = 2.0**(((temp_air(index_assi(inia))-tp_00)-25.0)/10.0)
1948	!
1949	kt(index_assi(inia)) = 0.7 * x_1 * 1.e6
1950	rt(index_assi(inia)) = 0.8 * x_1 / &
1951	( un + EXP(1.3*(temp_air(index_assi(inia))-tmax(index_assi(inia),jv))) )
1952	!
1953	vc(index_assi(inia)) = vcmax(index_assi(inia),jv) &
1954	* 0.39 * x_1 * water_lim(index_assi(inia)) / &
1955	! * 0.39 * x_1 / &
1956	( (1.0+EXP(0.3*(tmin(index_assi(inia),jv)-temp_air(index_assi(inia))))) &
1957	* (1.0+EXP(0.3*(temp_air(index_assi(inia))-topt(index_assi(inia),jv)))) )
1958	! @endcodeinc
1959	!
1960	ENDDO
1961	ENDIF
1962	!
1963	IF (nina .GT. 0) then
1964	! For the points where assimilation is not calculated, variables are initialized to specific values.
1965	DO inina=1,nina
1966	kt(index_non_assi(inina)) = zero
1967	rt(index_non_assi(inina)) = zero
1968	vc(index_non_assi(inina)) = zero
1969	ENDDO
1970	ENDIF
1971	!
1972	ELSE
1973	!
1974	! @addtogroup Photosynthesis
1975	! @{
1976	!
1977	! 2.3 Calculates temperature dependent parameters for C3 plants.\n
1978	!! \latexonly
1979	!! \input{diffuco_trans_co2_2.3.tex}
1980	!! \endlatexonly
1981	! @}
1982	!
1983	IF (nia .GT. 0) then
1984	!
1985	DO inia=1,nia
1986	!
1987	! @codeinc
1988	temp_lim(index_assi(inia)) = &
1989	(temp_air(index_assi(inia))-tmin(index_assi(inia),jv)) * &
1990	(temp_air(index_assi(inia))-tmax(index_assi(inia),jv))
1991	temp_lim(index_assi(inia)) = temp_lim(index_assi(inia)) / &
1992	(temp_lim(index_assi(inia))-(temp_air(index_assi(inia))-&
1993	topt(index_assi(inia),jv))**2)
1994	!
1995	Kc(index_assi(inia)) = 39.09 * EXP(.085*(temp_air(index_assi(inia))-tp_00))
1996	!
1997	Ko(index_assi(inia)) = 2.412 * 210000. &
1998	* EXP(.085*(temp_air(index_assi(inia))-tmin(index_assi(inia),jv))) / &
1999	(temp_air(index_assi(inia))-tmin(index_assi(inia),jv))
2000	!
2001	CP(index_assi(inia)) = 42. * EXP( 9.46*(temp_air(index_assi(inia))-tp_00-25.)/&
2002	temp_air(index_assi(inia)) )
2003	!
2004	vc(index_assi(inia)) = vcmax(index_assi(inia),jv) * &
2005	temp_lim(index_assi(inia)) * water_lim(index_assi(inia))
2006	! temp_lim(index_assi(inia))
2007	vj(index_assi(inia)) = vjmax(index_assi(inia),jv) * &
2008	temp_lim(index_assi(inia)) * water_lim(index_assi(inia))
2009	! temp_lim(index_assi(inia))
2010	! @endcodeinc
2011	!
2012	ENDDO
2013	!
2014	ENDIF
2015	!
2016	IF (nina .GT. 0) then
2017	! For the points where assimilation is not calculated, variables are initialized to specific values.
2018	DO inina=1,nina
2019	temp_lim(index_non_assi(inina)) = zero
2020	Kc(index_non_assi(inina)) = zero
2021	Ko(index_non_assi(inina)) = zero
2022	CP(index_non_assi(inina)) = zero
2023	!
2024	vc(index_non_assi(inina)) = zero
2025	vj(index_non_assi(inina)) = zero
2026	ENDDO
2027	ENDIF
2028	ENDIF ! C3/C4
2029	!
2030	! @addtogroup Photosynthesis
2031	! @{
2032	!
2033	! 2.4 Loop over LAI discretized levels to estimate assimilation and conductance\n
2034	! @}
2035	!
2036	!! The calculate(kjpindex) array is of type logical to indicate wether we have to sum over this LAI fixed level (the LAI of
2037	!! 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
2038	!! corresponding point is indexed in the index_calc array.
2039	DO jl = 1, nlai
2040	!
2041	nic=0
2042	calculate(:) = .FALSE.
2043	!
2044	IF (nia .GT. 0) then
2045	DO inia=1,nia
2046	calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) )
2047	IF ( calculate(index_assi(inia)) ) THEN
2048	nic=nic+1
2049	index_calc(nic)=index_assi(inia)
2050	ENDIF
2051	ENDDO
2052	ENDIF
2053	!
2054	! @addtogroup Photosynthesis
2055	! @{
2056	!
2057	! 2.4.1 Vmax is scaled into the canopy due to reduction of nitrogen
2058	!! (Johnson and Thornley,1984).\n
2059	!! \latexonly
2060	!! \input{diffuco_trans_co2_2.4.1.tex}
2061	!! \endlatexonly
2062	! @}
2063	!
2064	x_1 = ( un - .7_r_std * ( un - light(jv,jl) ) )
2065	!
2066	IF ( nic .GT. 0 ) THEN
2067	DO inic=1,nic
2068	! @codeinc
2069	vc2(index_calc(inic)) = vc(index_calc(inic)) * x_1
2070	vj2(index_calc(inic)) = vj(index_calc(inic)) * x_1
2071	! @endcodeinc
2072	ENDDO
2073	ENDIF
2074	!
2075	IF ( is_c4(jv) ) THEN
2076	!
2077	! @addtogroup Photosynthesis
2078	! @{
2079	!
2080	! 2.4.2 Assimilation for C4 plants (Collatz et al., 1992)\n
2081	!! \latexonly
2082	!! \input{diffuco_trans_co2_2.4.2.tex}
2083	!! \endlatexonly
2084	! @}
2085	!
2086	DO ji = 1, kjpindex
2087	assimi(ji) = zero
2088	ENDDO
2089	!
2090	IF (nic .GT. 0) THEN
2091	DO inic=1,nic
2092	! @codeinc
2093	x_1 = - ( vc2(index_calc(inic)) + 0.092 * 2.3* swdown(index_calc(inic)) * &
2094	ext_coef(jv) * light(jv,jl) )
2095	x_2 = vc2(index_calc(inic)) * 0.092 * 2.3 * swdown(index_calc(inic)) * &
2096	ext_coef(jv) * light(jv,jl)
2097	x_3 = ( -x_1 - sqrt( x_1x_1 - 4.0 xc4_1 * x_2 ) ) / (2.0*xc4_1)
2098	x_4 = - ( x_3 + kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * &
2099	1.0e-6 )
2100	x_5 = x_3 * kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * 1.0e-6
2101	assimi(index_calc(inic)) = ( -x_4 - sqrt( x_4x_4 - 4. xc4_2 * x_5 ) ) / (2.*xc4_2)
2102	assimi(index_calc(inic)) = assimi(index_calc(inic)) - &
2103	rt(index_calc(inic))
2104	! @endcodeinc
2105	ENDDO
2106	ENDIF
2107	ELSE
2108	!
2109	! @addtogroup Photosynthesis
2110	! @{
2111	!
2112	! 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
2113	!! \latexonly
2114	!! \input{diffuco_trans_co2_2.4.3.tex}
2115	!! \endlatexonly
2116	! @}
2117	!
2118	DO ji = 1, kjpindex
2119	assimi(ji) = zero
2120	ENDDO
2121	!
2122	IF (nic .GT. 0) THEN
2123	DO inic=1,nic
2124	! @codeinc
2125	x_1 = vc2(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) / &
2126	( leaf_ci(index_calc(inic),jv,jl) + Kc(index_calc(inic)) * &
2127	( un + 210000._r_std / Ko(index_calc(inic)) ) )
2128
2129	! The value 0.8855 is equal to 0.385*4.6/2.
2130	! The value 0.385 comes from MacMutrie et al. (1990), to calculate J based on the equation:
2131	! thetaJ^2 - (0.385theta+Jmax)J+O.385theta*Jmax=0 equation that comes
2132	! from Farqhuar and Wong (1984).
2133	! For C4 plants PAR is calculated from SWdown by dividing by 2 and multiplying by
2134	! by 4.6 to convert from W/m^2 to micromol/s/m^2.
2135	x_2 = .8855_r_stdswdown(index_calc(inic))ext_coef(jv)*light(jv,jl)
2136
2137	x_3 = x_2+vj2(index_calc(inic))
2138	x_4 = ( x_3 - sqrt( x_3x_3 - (4._r_std.7_r_stdx_2vj2(index_calc(inic))) ) ) / &
2139	(2._r_std*.7_r_std)
2140	x_5 = x_4 * leaf_ci(index_calc(inic),jv,jl) / &
2141	( 4.5_r_std * leaf_ci(index_calc(inic),jv,jl) + &
2142	10.5_r_std*CP(index_calc(inic)) )
2143	x_6 = MIN( x_1, x_5 )
2144	assimi(index_calc(inic)) = x_6 * ( un - CP(index_calc(inic))/&
2145	leaf_ci(index_calc(inic),jv,jl) ) - .011_r_std * vc2(index_calc(inic))
2146	! @endcodeinc
2147	ENDDO
2148	ENDIF
2149	ENDIF
2150	!
2151	IF (nic .GT. 0) THEN
2152	!
2153	DO inic=1,nic
2154	!
2155	! @addtogroup Photosynthesis
2156	! @{
2157	!
2158	!! 2.4.4 Estimatation of the stomatal conductance (Ball et al., 1987).\n
2159	!! \latexonly
2160	!! \input{diffuco_trans_co2_2.4.4.tex}
2161	!! \endlatexonly
2162	! @}
2163	!
2164	icinic=index_calc(inic)
2165	! gs(icinic) = water_lim(icinic) * &
2166	gs(icinic) = &
2167	( gsslope(jv) * assimi(icinic) * &
2168	air_relhum(icinic) / ccanopy(icinic) ) &
2169	+ gsoffset(jv)
2170	gs(icinic) = MAX( gs(icinic), gsoffset(jv) )
2171	ENDDO
2172	!
2173	DO inic=1,nic
2174	icinic=index_calc(inic)
2175	!
2176	! the new ci is calculated with
2177	! dci/dt=(ccanopy-ci)*gs/1.6-A
2178	! ci=ci+((ccanopy(icinic)-ci)*gs/1.6-&
2179	! assimi(icinic))*dtradia
2180	! we verify that ci is not out of possible values
2181	!
2182	ci_gs(icinic) = MIN( ccanopy(icinic), MAX( CP(icinic), &
2183	( ccanopy(icinic) - 1.6_r_std * assimi(icinic) / &
2184	gs(icinic) ) ) ) - leaf_ci(icinic,jv,jl)
2185	ENDDO
2186	DO inic=1,nic
2187	icinic=index_calc(inic)
2188	!to avoid some problem of numerical stability, the leaf_ci is bufferized
2189	!! [DISPENSABLE] This trick makes the code uncomprehensible:
2190	!! ci_gs = (Ca - 1.6*A)/gs -leaf_ci
2191	!! leaf_ci = leaf_ci + ci_gs/6. = leaf_ci + ((Ca - 1.6A)/gs - leaf_ci)/6. = leaf_ci5/6. + (Ca - 1.6A)/gs1/6.
2192	leaf_ci(icinic,jv,jl) = leaf_ci(icinic,jv,jl) + ci_gs(icinic)/6.
2193	ENDDO
2194	!
2195	DO inic=1,nic
2196	icinic=index_calc(inic)
2197	!
2198	! this might be the last level for which Ci is calculated. Store it for
2199	! initialization of the remaining levels of the Ci array.
2200	!
2201	leaf_ci_lowest(icinic) = leaf_ci(icinic,jv,jl)
2202	ENDDO
2203	!
2204	DO inic=1,nic
2205	icinic=index_calc(inic)
2206	!
2207	! @addtogroup Photosynthesis
2208	! @{
2209	!
2210	!! 2.4.5 Integration at the canopy level\n
2211	!! \latexonly
2212	!! \input{diffuco_trans_co2_2.4.5.tex}
2213	!! \endlatexonly
2214	! @}
2215	! total assimilation and conductance
2216	assimtot(icinic) = assimtot(icinic) + &
2217	assimi(icinic) * (laitab(jl+1)-laitab(jl))
2218	gstot(icinic) = gstot(icinic) + &
2219	gs(icinic) * (laitab(jl+1)-laitab(jl))
2220	!
2221	ilai(icinic) = jl
2222	!
2223	ENDDO
2224	!
2225	ENDIF
2226	!
2227	! keep stomatal conductance of topmost level
2228	!
2229	IF ( jl .EQ. 1 ) THEN
2230	!
2231	leaf_gs_top(:) = zero
2232	!
2233	IF ( nic .GT. 0 ) then
2234	DO inic=1,nic
2235	leaf_gs_top(index_calc(inic)) = gs(index_calc(inic))
2236	ENDDO
2237	ENDIF
2238	!
2239	ENDIF
2240	!
2241	IF (nia .GT. 0) THEN
2242	!
2243	DO inia=1,nia
2244	!
2245	IF ( .NOT. calculate(index_assi(inia)) ) THEN
2246	!
2247	! a) for plants that are doing photosynthesis, but whose LAI is lower
2248	! than the present LAI step, initialize it to the Ci of the lowest
2249	! canopy level
2250	!
2251	leaf_ci(index_assi(inia),jv,jl) = leaf_ci_lowest(index_assi(inia))
2252	!
2253	ENDIF
2254	!
2255	ENDDO
2256	!
2257	ENDIF
2258	!
2259	ENDDO ! loop over LAI steps
2260	!
2261	!! 2.5 Calculate resistances
2262	!
2263	IF (nia .GT. 0) THEN
2264	!
2265	DO inia=1,nia
2266	!
2267	iainia=index_assi(inia)
2268	!
2269	! conversion from mmol/m^2/s to m/s
2270	!
2271	gstot(iainia) = .0244(temp_air(iainia)/tp_00)&
2272	(1013./pb(iainia))*gstot(iainia)
2273	gstop(iainia) = .0244 * (temp_air(iainia)/tp_00)*&
2274	(1013./pb(iainia))leaf_gs_top(iainia)&
2275	laitab(ilai(iainia)+1)
2276	!
2277	rveget(iainia,jv) = un/gstop(iainia)
2278	!
2279	ENDDO
2280	!
2281	DO inia=1,nia
2282	!
2283	iainia=index_assi(inia)
2284	!
2285	! rstruct is the difference between rtot (=1./gstot) and rveget
2286	!
2287	! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
2288	!rstruct(iainia,jv) = un/gstot(iainia) - &
2289	! rveget(iainia,jv)
2290	rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
2291	rveget(iainia,jv), min_sechiba)
2292	!
2293	ENDDO
2294	!
2295	DO inia=1,nia
2296	!
2297	iainia=index_assi(inia)
2298	!
2299	!! wind is a global variable of the diffuco module.
2300	speed = MAX(min_wind, wind(index_assi(inia)))
2301	!
2302	! beta for transpiration
2303	!
2304	! Corrections Nathalie - 28 March 2006 - on advices of Fred Hourdin
2305	!! Introduction of a potentiometer rveg_pft to settle the rveg+rstruct sum problem in the coupled mode.
2306	!! rveg_pft=1 in the offline mode. rveg_pft is a global variable declared in the diffuco module.
2307	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2308	! (un - zqsvegrap(iainia)) * &
2309	! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
2310	! rstruct(iainia,jv))))
2311	!! Global resistance of the canopy to evaporation
2312	cresist=(un / (un + speed * q_cdrag(iainia) * &
2313	(rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
2314
2315
2316	vbeta3(iainia,jv) = veget_max(iainia,jv) * &
2317	(un - zqsvegrap(iainia)) * cresist + &
2318	!!$ ! Addition Nathalie - June 2006
2319	!!$ vbeta3(iainia,jv) = vbeta3(iainia,jv) + &
2320	! Corrections Nathalie - 09 November 2009 : veget => veget_max
2321	! MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
2322	MIN( vbeta23(iainia,jv), veget_max(iainia,jv) * &
2323	! zqsvegrap(iainia) * humrel(iainia,jv) * &
2324	zqsvegrap(iainia) * cresist )
2325	! Fin ajout Nathalie
2326	!
2327	! beta for assimilation. The difference is that surface
2328	! covered by rain (un - zqsvegrap(iainia)) is not taken into account
2329	! The factor 1.6 is conversion for H2O to CO2 conductance.
2330	! vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
2331	! (un / (un + q_cdrag(iainia) * &
2332	! (rveget(iainia,jv))))/1.6
2333	!
2334
2335	!! Global resistance of the canopy to CO2 transfer
2336	vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
2337	(un / (un + speed * q_cdrag(iainia) * &
2338	(rveget(iainia,jv) + rstruct(iainia,jv)))) / 1.6
2339	!
2340	! cimean is the "mean ci" calculated in such a way that assimilation
2341	! calculated in enerbil is equivalent to assimtot
2342	!
2343	cimean(iainia,jv) = ccanopy(iainia) - &
2344	assimtot(iainia) / &
2345	( vbetaco2(iainia,jv)/veget_max(iainia,jv) * &
2346	rau(iainia) * speed * q_cdrag(iainia))
2347	!
2348	ENDDO
2349	!
2350	ENDIF
2351	!
2352	END DO ! loop over vegetation types
2353	!
2354	IF (long_print) WRITE (numout,*) ' diffuco_trans_co2 done '
2355
2356
2357	END SUBROUTINE diffuco_trans_co2
2358
2359
2360	!! ================================================================================================================================
2361	!! SUBROUTINE : diffuco_comb
2362	!!
2363	!>\BRIEF This routine combines the previous partial beta
2364	!! coefficients and calculates the total alpha and complete beta coefficients.
2365	!!
2366	!! DESCRIPTION : Those integrated coefficients are used to calculate (in enerbil.f90) the total evapotranspiration
2367	!! from the grid-cell. \n
2368	!!
2369	!! In the case that air is more humid than surface, dew deposition can occur (negative latent heat flux).
2370	!! In this instance, for temperature above zero, all of the beta coefficients are set to 0, except for
2371	!! interception (vbeta2) and bare soil (vbeta4 with zero soil resistance). The amount of water that is
2372	!! intercepted by leaves is calculated based on the value of LAI of the surface. In the case of freezing
2373	!! temperatures, water is added to the snow reservoir, and so vbeta4 and vbeta2 are set to 0, and the
2374	!! total vbeta and valpha is set to 1.\n
2375	!!
2376	!! \latexonly
2377	!! \input{diffucocomb1.tex}
2378	!! \endlatexonly
2379	!!
2380	!! The beta and alpha coefficients are initially set to 1.
2381	!! \latexonly
2382	!! \input{diffucocomb2.tex}
2383	!! \endlatexonly
2384	!!
2385	!! If snow is lower than the critical value:
2386	!! \latexonly
2387	!! \input{diffucocomb3.tex}
2388	!! \endlatexonly
2389	!! If in the presence of dew:
2390	!! \latexonly
2391	!! \input{diffucocomb4.tex}
2392	!! \endlatexonly
2393	!!
2394	!! Determine where the water goes (soil, vegetation, or snow)
2395	!! when air moisture exceeds saturation.
2396	!! \latexonly
2397	!! \input{diffucocomb5.tex}
2398	!! \endlatexonly
2399	!!
2400	!! If it is not freezing dew is put into the interception reservoir and onto the bare soil. If it is freezing,
2401	!! water is put into the snow reservoir.
2402	!! Now modify valpha and vbetas where necessary: for soil and snow
2403	!! \latexonly
2404	!! \input{diffucocomb6.tex}
2405	!! \endlatexonly
2406	!!
2407	!! and for vegetation
2408	!! \latexonly
2409	!! \input{diffucocomb7.tex}
2410	!! \endlatexonly
2411	!!
2412	!! Then compute part of dew that can be intercepted by leafs.
2413	!!
2414	!! There will be no transpiration when air moisture is too high, under any circumstance
2415	!! \latexonly
2416	!! \input{diffucocomb8.tex}
2417	!! \endlatexonly
2418	!!
2419	!! There will also be no interception loss on bare soil, under any circumstance.
2420	!! \latexonly
2421	!! \input{diffucocomb9.tex}
2422	!! \endlatexonly
2423	!!
2424	!! The flowchart details the 'decision tree' which underlies the module.
2425	!!
2426	!! RECENT CHANGE(S): None
2427	!!
2428	!! MAIN OUTPUT VARIABLE(S): vbeta1, vbeta4, humrel, vbeta2, vbeta3, valpha, vbeta
2429	!!
2430	!! REFERENCE(S) :
2431	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
2432	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
2433	!! Circulation Model. Journal of Climate, 6, pp.248-273
2434	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influences de l'irrigation
2435	!! sur le cycle de l'eau, PhD Thesis, available from:
2436	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
2437	!!
2438	!! FLOWCHART :
2439	!! \latexonly
2440	!! \includegraphics[scale=0.25]{diffuco_comb_flowchart.png}
2441	!! \endlatexonly
2442	!! \n
2443	!_ ================================================================================================================================
2444
2445	SUBROUTINE diffuco_comb (kjpindex, dtradia, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
2446	& snow, veget, lai, vbeta1, vbeta2, vbeta3 , vbeta4, valpha, vbeta, qsintmax)
2447
2448	! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006
2449
2450	!! 0. Variable and parameter declaration
2451
2452	!! 0.1 Input variables
2453
2454	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
2455	REAL(r_std), INTENT (in) :: dtradia !! Time step (s)
2456	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Air Density (kg m^{-3})
2457	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind speed (m s^{-1})
2458	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Nortward Lowest level wind speed (m s^{-1})
2459	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Product of Surface drag coefficient and wind speed (m s^{-1})
2460	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure (Pa)
2461	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific air humidity (kg kg^{-1})
2462	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature (K)
2463	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lower air temperature (K)
2464	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass (kg)
2465	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type (fraction)
2466	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
2467	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation (kg m^{-2})
2468
2469	!! 0.2 Output variables
2470
2471	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: valpha !! Total Alpha coefficient (-)
2472	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient (-)
2473
2474	!! 0.3 Modified variables
2475
2476	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation (-)
2477	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation (-)
2478	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress (within range 0 to 1)
2479	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for interception loss (-)
2480	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration (-)
2481
2482	!! 0.4 Local variables
2483
2484	INTEGER(i_std) :: ji, jv
2485	REAL(r_std) :: zevtest, zsoil_moist, zrapp
2486	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum
2487	REAL(r_std), DIMENSION(kjpindex) :: vegetsum, vegetsum2
2488	REAL(r_std), DIMENSION(kjpindex) :: qsatt
2489	LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow
2490	REAL(r_std) :: coeff_dew_veg
2491	!_ ================================================================================================================================
2492
2493	!! \latexonly
2494	!! \input{diffucocomb1.tex}
2495	!! \endlatexonly
2496	vbeta2sum(:) = zero
2497	vbeta3sum(:) = zero
2498	DO jv = 1, nvm
2499	vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv)
2500	vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv)
2501	ENDDO
2502
2503	!! 1. The beta and alpha coefficients are initially set to 1.
2504
2505	!! \latexonly
2506	!! \input{diffucocomb2.tex}
2507	!! \endlatexonly
2508	vbeta(:) = un
2509	valpha(:) = un
2510
2511
2512	!! 2. if snow is lower than the critical value:
2513
2514	!! \latexonly
2515	!! \input{diffucocomb3.tex}
2516	!! \endlatexonly
2517	DO ji = 1, kjpindex
2518
2519	IF (snow(ji) .LT. snowcri) THEN
2520
2521	vbeta(ji) = vbeta4(ji) + vbeta2sum(ji) + vbeta3sum(ji)
2522
2523	IF (vbeta(ji) .LT. min_sechiba) THEN
2524	vbeta(ji) = zero
2525	END IF
2526
2527	END IF
2528
2529	ENDDO
2530
2531
2532	!! 3 If we are in presence of dew:
2533
2534	! for vectorization: some arrays
2535
2536	!! \latexonly
2537	!! \input{diffucocomb4.tex}
2538	!! \endlatexonly
2539	vegetsum(:) = zero
2540
2541	DO jv = 1, nvm
2542
2543	vegetsum(:) = vegetsum(:) + veget(:,jv)
2544
2545	ENDDO
2546
2547	vegetsum2(:) = zero
2548
2549	DO jv = 2, nvm
2550
2551	vegetsum2(:) = vegetsum2(:) + veget(:,jv)
2552
2553	ENDDO
2554
2555	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
2556
2557
2558	!! 9.3.1 Determine where the water goes
2559	!! Determine where the water goes (soil, vegetation, or snow)
2560	!! when air moisture exceeds saturation.
2561	!! \latexonly
2562	!! \input{diffucocomb5.tex}
2563	!! \endlatexonly
2564	toveg(:) = .FALSE.
2565	tosnow(:) = .FALSE.
2566	DO ji = 1, kjpindex
2567	IF ( qsatt(ji) .LT. qair(ji) ) THEN
2568	IF (temp_air(ji) .GT. tp_00) THEN
2569
2570	!! 9.3.1.1 If it is not freezing,
2571	!! If it is not freezing dew is put into the
2572	!! interception reservoir and onto the bare soil.
2573	toveg(ji) = .TRUE.
2574	ELSE
2575
2576	!! 9.3.1.2 If it is freezing,
2577	!! If it is freezing water is put into the
2578	!! snow reservoir.
2579	tosnow(ji) = .TRUE.
2580	ENDIF
2581	ENDIF
2582	END DO
2583
2584
2585	!! 9.3.1.3 Now modify valpha and vbetas where necessary.
2586
2587	!! 9.3.1.3.1 Soil and snow (2d)
2588	!! \latexonly
2589	!! \input{diffucocomb6.tex}
2590	!! \endlatexonly
2591	DO ji = 1, kjpindex
2592	IF ( toveg(ji) ) THEN
2593	vbeta1(ji) = zero
2594	vbeta4(ji) = veget(ji,1)
2595	! Correction Nathalie - le 13-03-2006: le vbeta ne sera calcule qu'une fois tous les vbeta2 redefinis
2596	!vbeta(ji) = vegetsum(ji)
2597	vbeta(ji) = vbeta4(ji)
2598	valpha(ji) = un
2599	ENDIF
2600	IF ( tosnow(ji) ) THEN
2601	vbeta1(ji) = un
2602	vbeta4(ji) = zero
2603	vbeta(ji) = un
2604	valpha(ji) = un
2605	ENDIF
2606	ENDDO
2607
2608	!! 9.3.1.3.2 Vegetation (3d)
2609	!! \latexonly
2610	!! \input{diffucocomb7.tex}
2611	!! \endlatexonly
2612	DO jv = 1, nvm
2613
2614	DO ji = 1, kjpindex
2615
2616	! Correction Nathalie - 13-03-2006 / si qsintmax=0, vbeta2=0
2617	IF ( toveg(ji) ) THEN
2618	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
2619
2620	! Compute part of dew that can be intercepted by leafs.
2621	IF ( lai(ji,jv) .GT. min_sechiba) THEN
2622	IF (lai(ji,jv) .GT. 1.5) THEN
2623	coeff_dew_veg= &
2624	& 0.000017lai(ji,jv)*5 &
2625	& - 0.000824lai(ji,jv)*4 &
2626	& + 0.014843lai(ji,jv)*3 &
2627	& - 0.110112lai(ji,jv)*2 &
2628	& + 0.205673*lai(ji,jv) &
2629	& + 0.887773
2630	ELSE
2631	coeff_dew_veg=un
2632	ENDIF
2633	ELSE
2634	coeff_dew_veg=zero
2635	ENDIF
2636	vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv)
2637	!vbeta2(ji,jv) = veget(ji,jv)
2638	ELSE
2639	vbeta2(ji,jv) = zero
2640	ENDIF
2641	vbeta(ji) = vbeta(ji) + vbeta2(ji,jv)
2642	ENDIF
2643	IF ( tosnow(ji) ) vbeta2(ji,jv) = zero
2644
2645	ENDDO
2646
2647	ENDDO
2648
2649	!! 9.3.2a There will be no transpiration when air moisture is too high, under any circumstance
2650	!! \latexonly
2651	!! \input{diffucocomb8.tex}
2652	!! \endlatexonly
2653	DO jv = 1, nvm
2654	DO ji = 1, kjpindex
2655	IF ( qsatt(ji) .LT. qair(ji) ) THEN
2656	vbeta3(ji,jv) = zero
2657	humrel(ji,jv) = zero
2658	ENDIF
2659	ENDDO
2660	ENDDO
2661
2662
2663	!! 9.3.2b There will also be no interception loss on bare soil, under any circumstance.
2664	!! \latexonly
2665	!! \input{diffucocomb9.tex}
2666	!! \endlatexonly
2667	DO ji = 1, kjpindex
2668	IF ( qsatt(ji) .LT. qair(ji) ) THEN
2669	vbeta2(ji,1) = zero
2670	ENDIF
2671	ENDDO
2672
2673	IF (long_print) WRITE (numout,*) ' diffuco_comb done '
2674
2675	END SUBROUTINE diffuco_comb
2676
2677
2678	!! ================================================================================================================================
2679	!! SUBROUTINE : diffuco_raerod
2680	!!
2681	!>\BRIEF Computes the aerodynamic resistance, for cases in which the
2682	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ and when the flag
2683	!! 'ldq_cdrag_from_gcm' is set to TRUE
2684	!!
2685	!! DESCRIPTION : Simply computes the aerodynamic resistance, for cases in which the
2686	!! surface drag coefficient is provided by the coupled atmospheric model LMDZ. If the surface drag coefficient
2687	!! is not provided by the LMDZ or signalled by the flag 'ldq_cdrag_from_gcm' set to FALSE, then the subroutine
2688	!! diffuco_aero is called instead of this one.
2689	!!
2690	!! Calculation of the aerodynamic resistance, for diganostic purposes. First calculate wind speed:
2691	!! \latexonly
2692	!! \input{diffucoaerod1.tex}
2693	!! \endlatexonly
2694	!!
2695	!! next calculate ::raero
2696	!! \latexonly
2697	!! \input{diffucoaerod2.tex}
2698	!! \endlatexonly
2699	!!
2700	!! RECENT CHANGE(S): None
2701	!!
2702	!! MAIN OUTPUT VARIABLE(S): ::raero
2703	!!
2704	!! REFERENCE(S) :
2705	!! - de Noblet-DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
2706	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
2707	!! Circulation Model. Journal of Climate, 6, pp.248-273
2708	!! - Guimberteau, M, 2010. ModÃ©lisation de l'hydrologie continentale et influence de l'irrigation
2709	!! sur le cycle de l'eau, PhD Thesis, available from:
2710	!! http://www.sisyphe.upmc.fr/~guimberteau/docs/manuscrit_these.pdf
2711	!!
2712	!! FLOWCHART : None
2713	!! \n
2714	!_ ================================================================================================================================
2715
2716	SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
2717
2718	IMPLICIT NONE
2719
2720	!! 0. Variable and parameter declaration
2721
2722	!! 0.1 Input variables
2723
2724	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (-)
2725	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Eastward Lowest level wind velocity (m s^{-1})
2726	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Northward Lowest level wind velocity (m s^{-1})
2727	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Product of Surface drag coefficient and wind speed (m s^{-1})
2728
2729	!! 0.2 Output variables
2730
2731	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance (s m^{-1})
2732
2733	!! 0.3 Modified variables
2734
2735	!! 0.4 Local variables
2736
2737	INTEGER(i_std) :: ji !! (-)
2738	REAL(r_std) :: speed !! (m s^{-1})
2739	!_ ================================================================================================================================
2740
2741	!! 1. Simple calculation of the aerodynamic resistance, for diganostic purposes.
2742
2743	DO ji=1,kjpindex
2744
2745	!! \latexonly
2746	!! \input{diffucoaerod1.tex}
2747	!! \endlatexonly
2748	speed = MAX(min_wind, wind(ji))
2749
2750	!! \latexonly
2751	!! \input{diffucoaerod2.tex}
2752	!! \endlatexonly
2753	raero(ji) = un / (q_cdrag(ji)*speed)
2754
2755	ENDDO
2756
2757	END SUBROUTINE diffuco_raerod
2758
2759
2760	END MODULE diffuco

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

Context Navigation

source: perso/abdelouhab.djerrah/ORCHIDEE/src_sechiba/diffuco.f90 @ 854

Download in other formats: