Context Navigation

source: tags/ORCHIDEE/src_sechiba/diffuco.f90 @ 6

Last change on this file since 6 was 6, checked in by orchidee, 14 years ago
import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 71.2 KB

Line
1	!!
2	!! This module computes diffusion coefficients for continental points.
3	!!
4	!! @author Marie-Alice Foujols and Jan Polcher
5	!! @Version : $Revision: 1.35 $, $Date: 2010/04/07 09:16:40 $
6	!!
7	!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/diffuco.f90,v 1.35 2010/04/07 09:16:40 ssipsl Exp $
8	!! IPSL (2006)
9	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10	!!
11	MODULE diffuco
12
13	! modules used :
14	USE constantes
15	USE constantes_veg
16	USE sechiba_io
17	USE ioipsl
18	USE constantes_co2
19	USE parallel
20	! USE WRITE_FIELD_p
21
22	IMPLICIT NONE
23
24	! public routines :
25	! diffuco_main only
26	PRIVATE
27	PUBLIC :: diffuco_main,diffuco_clear
28
29	!
30	! variables used inside diffuco module : declaration and initialisation
31	!
32	INTEGER(i_std), PARAMETER :: nlai = 20
33	LOGICAL, SAVE :: l_first_diffuco = .TRUE. !! Initialisation has to be done one time
34	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
35	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: leaf_ci !! intercellular CO2 concentration (ppm)
36	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rstruct !! architectural resistance
37	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: raero !! Aerodynamic resistance
38	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsatt !! Surface saturated humidity
39	!! Nathalie le 28 mars 2006 - sur proposition de Fred Hourdin, ajout
40	!! d'un potentiometre pour regler la resistance de la vegetation
41	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rveg_pft
42	! MM
43	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: wind !! Wind norm
44
45	CONTAINS
46
47	!!
48	!! Main routine for diffuco module.
49	!! - called only one time for initialisation
50	!! - called every time step
51	!! - called one more time at last time step for writing _restart_ file
52	!!
53	!! Algorithm:
54	!! - call diffuco_aero for aerodynamic transfer coeficient
55	!! - call diffuco_snow for partial beta coefficient : sublimation
56	!! - call diffuco_inter for partial beta coefficient : interception for each type of vegetation
57	!! - call diffuco_bare for partial beta coefficient : bare soil
58	!! - call diffuco_trans for partial beta coefficient : transpiration for each type of vegetation
59	!! - call diffuco_comb for alpha and beta coefficient
60	!!
61	!! @call diffuco_aero
62	!! @call diffuco_snow
63	!! @call diffuco_inter
64	!! @call diffuco_bare
65	!! @call diffuco_trans
66	!! @call diffuco_comb
67	!!
68	SUBROUTINE diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
69	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
70	! & zlev, z0, roughheight, temp_sol, temp_air, rau, q_cdrag, qsurf, qair, pb, &
71	& zlev, z0, roughheight, temp_sol, temp_air, rau, q_cdrag, qsurf, qair, q2m, t2m, pb, &
72	& rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
73	& swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
74	& vbeta , valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
75
76	! interface description:
77	! input scalar
78	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
79	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
80	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
81	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
82	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier
83	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
84	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for restart file to read
85	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for restart file to write
86	! input fields
87	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
88	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map
89	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
90	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
91	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
92	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0 !! Surface roughness
93	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective height for roughness
94	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature in Kelvin
95	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature in Kelvin
96	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
97	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific humidity
98	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
99	! Ajout Nathalie - declaration q2m
100	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity
101	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature
102	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass
103	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Surface level pressure
104	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rsol !! Bare soil evaporation resistance
105	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! Beta factor for bare soil evaporation
106	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation
107	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,...
108	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,...
109	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+cities+...
110	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
111	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
112	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration inside the canopy
113	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
114	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI->infty)
115	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index
116	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
117	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
118	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in):: assim_param !! min+max+opt temps, vcmax, vjmax for photosynthesis
119	! modified fields
120	REAL(r_std),DIMENSION (kjpindex, nvm), INTENT (inout) :: humrel !! Moisture stress
121	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag ! Aerodynamic conductance
122	! output fields
123	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient
124	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: valpha !! Total alpha coefficient
125	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation
126	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation
127	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE: Beta for CO2
128	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss
129	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for transpiration
130	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Surface resistance for the vegetatuon
131	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE: mean intercellular ci (see enerbil)
132	! Local
133	! AJout Nathalie - Juin 2006
134	!! Beta for fraction of wetted foliage that will transpire
135	REAL(r_std),DIMENSION (kjpindex,nvm) :: vbeta23
136
137	INTEGER(i_std) :: ilai
138	CHARACTER(LEN=4) :: laistring
139
140	! do initialisation if needed
141
142	IF (l_first_diffuco) THEN
143
144	!Config Key = CDRAG_FROM_GCM
145	!Config Desc = Keep cdrag coefficient from gcm.
146	!Config Def = TRUE if q_cdrag on initialization is non zero
147	!Config Help = Set to .TRUE. if you want q_cdrag coming from GCM.
148	!Congig Keep cdrag coefficient from gcm for latent and sensible heat fluxes.
149	IF ( ABS(MAXVAL(q_cdrag)) .LE. EPSILON(q_cdrag)) THEN
150	ldq_cdrag_from_gcm = .FALSE.
151	ELSE
152	ldq_cdrag_from_gcm = .TRUE.
153	ENDIF
154	!MM q_cdrag is always 0 on initialization ??
155	CALL getin_p('CDRAG_from_GCM', ldq_cdrag_from_gcm)
156
157	WRITE(numout,*) "ldq_cdrag_from_gcm = ",ldq_cdrag_from_gcm
158
159	IF (long_print) WRITE (numout,*) ' call diffuco_init '
160
161	! If cdrag is
162	CALL diffuco_init(kjit, ldrestart_read, kjpindex, index, rest_id, q_cdrag)
163
164	RETURN
165
166	ENDIF
167	!
168	! prepares restart file for the next simulation
169	!
170	IF (ldrestart_write) THEN
171
172	IF (long_print) WRITE (numout,*) ' we have to complete restart file with DIFFUCO variables '
173
174	var_name= 'rstruct'
175	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, rstruct, 'scatter', nbp_glo, index_g)
176
177	var_name= 'raero'
178	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, raero, 'scatter', nbp_glo, index_g)
179
180	var_name= 'qsatt'
181	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, qsatt, 'scatter', nbp_glo, index_g)
182
183	! the following variable is written only if CO2 was calculated
184	IF ( control%ok_co2 ) THEN
185
186	DO ilai = 1, nlai
187
188	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
189	write(laistring,'(i4)') ilai
190	laistring=ADJUSTL(laistring)
191	var_name='leaf_ci_'//laistring(1:LEN_TRIM(laistring))
192
193	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, leaf_ci(:,:,ilai), 'scatter', nbp_glo, index_g)
194
195	ENDDO
196
197	ENDIF
198
199	IF (.NOT.ldq_cdrag_from_gcm) THEN
200	var_name= 'cdrag'
201	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, q_cdrag, 'scatter', nbp_glo, index_g)
202	ENDIF
203
204	RETURN
205
206	END IF
207	! MM
208	wind(:) = SQRT (u(:)u(:) + v(:)v(:))
209
210	!
211	! calculs des differents coefficients
212	!
213	IF (.NOT.ldq_cdrag_from_gcm) THEN
214	CALL diffuco_aero (kjpindex, kjit, u, v, zlev, z0, roughheight, veget_max, temp_sol, temp_air, &
215	& qsurf, qair, q_cdrag)
216	ENDIF
217	CALL diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
218
219	!
220	! An estimation of the satturated humidity at the surface
221	!
222
223	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
224
225	!
226	! beta coefficient for sublimation
227	!
228
229	CALL diffuco_snow (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, &
230	& snow, frac_nobio, totfrac_nobio, snow_nobio, vbeta1)
231
232	!
233	! beta coefficient for interception
234	!
235
236	! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23
237	!CALL diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
238	! & qsintveg, qsintmax, rstruct, vbeta2)
239	CALL diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
240	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
241
242	!
243	! beta coefficient for bare soil
244	!
245
246	CALL diffuco_bare (kjpindex, dtradia, u, v, q_cdrag, rsol, evap_bare_lim, evapot, humrel, veget, vbeta4)
247
248	!
249	! beta coefficient for transpiration
250	!
251
252	IF ( control%ok_co2 ) THEN
253
254	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
255	! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23
256	!CALL diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
257	! assim_param, ccanopy, &
258	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
259
260	CALL diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, &
261	assim_param, ccanopy, &
262	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
263
264	ELSE
265
266	! Correction Nathalie - Juin 2006 - introduction d'un terme vbeta23
267	!CALL diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
268	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
269
270	CALL diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
271	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
272
273	ENDIF
274
275	!
276	! combination of coefficient : alpha and beta coefficient
277	!
278
279	! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006
280	CALL diffuco_comb (kjpindex, dtradia, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, &
281	& veget, lai, vbeta1, vbeta2, vbeta3, vbeta4, valpha, vbeta, qsintmax)
282
283	IF ( .NOT. almaoutput ) THEN
284	CALL histwrite(hist_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
285	CALL histwrite(hist_id, 'raero', kjit, raero, kjpindex, index)
286	! Ajouts Nathalie - novembre 2006
287	CALL histwrite(hist_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
288	CALL histwrite(hist_id, 'Wind', kjit, wind, kjpindex, index)
289	! Fin ajouts Nathalie
290	!MM
291	CALL histwrite(hist_id, 'qsatt', kjit, qsatt, kjpindex, index)
292
293	IF ( hist2_id > 0 ) THEN
294	CALL histwrite(hist2_id, 'rstruct', kjit, rstruct, kjpindex*nvm, indexveg)
295	CALL histwrite(hist2_id, 'raero', kjit, raero, kjpindex, index)
296	CALL histwrite(hist2_id, 'cdrag', kjit, q_cdrag, kjpindex, index)
297	CALL histwrite(hist2_id, 'Wind', kjit, wind, kjpindex, index)
298	CALL histwrite(hist2_id, 'qsatt', kjit, qsatt, kjpindex, index)
299	ENDIF
300	ELSE
301
302	ENDIF
303	!
304	!
305	IF (long_print) WRITE (numout,*) ' diffuco_main done '
306
307	END SUBROUTINE diffuco_main
308
309	!! Algorithm:
310	!! - dynamic allocation for local array
311	!!
312	SUBROUTINE diffuco_init(kjit, ldrestart_read, kjpindex, index, rest_id, q_cdrag)
313
314	! interface description
315	! input scalar
316	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
317	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read
318	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
319	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in):: index !! Indeces of the points on the map
320	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
321	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag
322	! input fields
323	! output scalar
324	! output fields
325
326	! local declaration
327	INTEGER(i_std) :: ier, jv
328	INTEGER(i_std) :: ilai
329	CHARACTER(LEN=4) :: laistring
330	REAL(r_std),DIMENSION (kjpindex) :: temp
331
332	!
333	! initialisation
334	!
335	IF (l_first_diffuco) THEN
336	l_first_diffuco=.FALSE.
337	ELSE
338	WRITE (numout,*) ' l_first_diffuco false . we stop '
339	STOP 'diffuco_init'
340	ENDIF
341
342	! allocate only if CO2 is calculated
343	IF ( control%ok_co2 ) THEN
344
345	ALLOCATE (leaf_ci(kjpindex,nvm,nlai),stat=ier)
346	IF (ier.NE.0) THEN
347	WRITE (numout,) ' error in leaf_ci allocation. We stop. We need kjpindexnvm*nlai words = ',&
348	kjpindexnvmnlai
349	STOP 'diffuco_init'
350	END IF
351
352	ENDIF
353
354	ALLOCATE (rstruct(kjpindex,nvm),stat=ier)
355	IF (ier.NE.0) THEN
356	WRITE (numout,*) ' error in rstruct allocation. We stop. We need kjpindex x nvm words = ' ,kjpindex,' x ' ,nvm,&
357	& ' = ',kjpindex*nvm
358	STOP 'diffuco_init'
359	END IF
360	ALLOCATE (raero(kjpindex),stat=ier)
361	IF (ier.NE.0) THEN
362	WRITE (numout,*) ' error in raero allocation. We stop. We need kjpindex x nvm words = ', kjpindex
363	STOP 'diffuco_init'
364	END IF
365
366	ALLOCATE (qsatt(kjpindex),stat=ier)
367	IF (ier.NE.0) THEN
368	WRITE (numout,*) ' error in qsatt allocation. We stop. We need kjpindex x nvm words = ', kjpindex
369	STOP 'diffuco_init'
370	END IF
371
372	ALLOCATE (wind(kjpindex),stat=ier)
373	IF (ier.NE.0) THEN
374	WRITE (numout,*) ' error in wind allocation. We stop. We need kjpindex x nvm words = ', kjpindex
375	STOP 'diffuco_init'
376	END IF
377
378	IF (ldrestart_read) THEN
379
380	IF (long_print) WRITE (numout,*) ' we have to read a restart file for DIFFUCO variables'
381
382	var_name='rstruct'
383	CALL ioconf_setatt('UNITS', 's/m')
384	CALL ioconf_setatt('LONG_NAME','Structural resistance')
385	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., rstruct, "gather", nbp_glo, index_g)
386	!
387	IF ( MINVAL(rstruct) .EQ. MAXVAL(rstruct) .AND. MAXVAL(rstruct) .EQ. val_exp ) THEN
388	!
389	DO jv = 1, nvm
390	rstruct(:,jv) = rstruct_const(jv)
391	ENDDO
392
393	ENDIF
394
395	var_name='raero' ;
396	CALL ioconf_setatt('UNITS', 's/m')
397	CALL ioconf_setatt('LONG_NAME','Aerodynamic resistance')
398	IF ( ok_var(var_name) ) THEN
399	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
400	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
401	raero(:) = temp(:)
402	ENDIF
403	ENDIF
404	!
405	var_name='qsatt' ;
406	CALL ioconf_setatt('UNITS', 'g/g')
407	CALL ioconf_setatt('LONG_NAME','Surface saturated humidity')
408	IF ( ok_var(var_name) ) THEN
409	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
410	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
411	qsatt(:) = temp(:)
412	ENDIF
413	ENDIF
414
415	! the following variable is read only if CO2 is calculated
416	IF ( control%ok_co2 ) THEN
417
418	CALL ioconf_setatt('UNITS', 'ppm')
419	CALL ioconf_setatt('LONG_NAME','Leaf CO2')
420
421	DO ilai = 1, nlai
422
423	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
424	write(laistring,'(i4)') ilai
425	laistring=ADJUSTL(laistring)
426	var_name='leaf_ci_'//laistring(1:LEN_TRIM(laistring))
427
428	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE.,leaf_ci(:,:,ilai), "gather", nbp_glo, index_g)
429
430	ENDDO
431	!
432	!Config Key = DIFFUCO_LEAFCI
433	!Config Desc = Initial leaf CO2 level if not found in restart
434	!Config Def = 233.
435	!Config Help = The initial value of leaf_ci if its value is not found
436	!Config in the restart file. This should only be used if the model is
437	!Config started without a restart file.
438	!
439	CALL setvar_p (leaf_ci, val_exp,'DIFFUCO_LEAFCI', 233._r_std)
440	ENDIF
441	!
442	IF (.NOT.ldq_cdrag_from_gcm) THEN
443	var_name= 'cdrag'
444	CALL ioconf_setatt('LONG_NAME','Drag coefficient for LE and SH')
445	CALL ioconf_setatt('UNITS', '-')
446	IF ( ok_var(var_name) ) THEN
447	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., temp, "gather", nbp_glo, index_g)
448	IF (MINVAL(temp) < MAXVAL(temp) .OR. MAXVAL(temp) .NE. val_exp) THEN
449	q_cdrag(:) = temp(:)
450	ENDIF
451	ENDIF
452
453	ENDIF
454
455	ENDIF
456
457	!
458	! Ajouts Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin
459	!
460	ALLOCATE (rveg_pft(nvm),stat=ier)
461	IF (ier.NE.0) THEN
462	WRITE (numout,*) &
463	' error in rveg_pft allocation. We stop. We need nvm words = ', nvm
464	STOP 'diffuco_init'
465	END IF
466	!
467	!Config Key = RVEG_PFT
468	!Config Desc = Artificial parameter to increase or decrease canopy resistance.
469	!Config Def = 1.
470	!Config Help = This parameter is set by PFT.
471	rveg_pft(:) = 1.
472	CALL getin_p('RVEG_PFT', rveg_pft)
473
474	WRITE(numout,*) 'DANS DIFFUCO_INIT , RVEG_PFT=',rveg_pft
475
476	IF (long_print) WRITE (numout,*) ' diffuco_init done '
477
478	END SUBROUTINE diffuco_init
479	SUBROUTINE diffuco_clear()
480
481	l_first_diffuco=.TRUE.
482
483	IF (ALLOCATED (leaf_ci)) DEALLOCATE (leaf_ci)
484	IF (ALLOCATED (rstruct)) DEALLOCATE (rstruct)
485	IF (ALLOCATED (raero)) DEALLOCATE (raero)
486	IF (ALLOCATED (rveg_pft)) DEALLOCATE (rveg_pft)
487
488	END SUBROUTINE diffuco_clear
489
490	!! This routine computes aerothermic coefficient if required
491	!! see logical ldq_cdrag_from_gcm
492	!!
493	SUBROUTINE diffuco_aero (kjpindex, kjit, u, v, zlev, z0, roughheight, veget_max, temp_sol, temp_air, &
494	& qsurf, qair, q_cdrag)
495
496	! interface description
497	! input scalar
498	INTEGER(i_std), INTENT(in) :: kjpindex, kjit !! Domain size
499	! input fields
500	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
501	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
502	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
503	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: z0 !! Surface roughness
504	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: roughheight !! Effective roughness height
505	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Fraction of vegetation type
506	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature in Kelvin
507	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Lowest level temperature in Kelvin
508	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsurf !! near surface specific humidity
509	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
510
511	! output scalar
512	! output fields
513	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: q_cdrag !! Surface drag ! Aerodynamic conductance
514
515	! local declaration
516	INTEGER(i_std) :: ji, jv
517	REAL(r_std) :: speed, zg, zdphi, ztvd, ztvs, zdu2
518	REAL(r_std) :: zri, cd_neut, zscf, cd_tmp
519
520	! initialisation
521
522	! test if we have to work with q_cdrag or to calcul it
523
524	DO ji=1,kjpindex
525	!
526	! 1. computes wind speed
527	!
528	speed = wind(ji)
529	!
530	! 2. computes geopotentiel
531	!
532	zg = zlev(ji) * cte_grav
533	zdphi = zg/cp_air
534	!
535	! 3. virtual air temperature at the surface
536	!
537	ztvd = (temp_air(ji) + zdphi / (un + rvtmp2 * qair(ji))) * (un + retv * qair(ji))
538	!
539	! 4. virtual surface temperature
540	!
541	ztvs = temp_sol(ji) * (un + retv * qsurf(ji))
542	!
543	! 5. squared wind shear
544	!
545	zdu2 = MAX(cepdu2,speed**2)
546	!
547	! 6. Richardson number
548	!
549	zri = zg * (ztvd - ztvs) / (zdu2 * ztvd)
550	zri = MAX(MIN(zri,5.),-5.)
551	!
552	! 7. Computing the drag coefficient
553	! We add the add the height of the vegetation to the level height to take into account
554	! that the level seen by the vegetation is actually the top of the vegetation. Then we
555	! we can subtract the displacement height.
556	!
557	cd_neut = (ct_karman / LOG( (zlev(ji) + roughheight(ji)) / z0(ji) )) ** 2
558	!
559	! 7.1 Stable case
560	!
561	IF (zri .GE. zero) THEN
562	zscf = SQRT(un + cd * ABS(zri))
563	cd_tmp=cd_neut/(un + trois * cb * zri * zscf)
564	ELSE
565	!
566	! 7.2 Unstable case
567	!
568	zscf = un / (un + trois * cb * cc * cd_neut * SQRT(ABS(zri) * &
569	& ((zlev(ji) + roughheight(ji)) / z0(ji))))
570	cd_tmp=cd_neut * (un - trois * cb * zri * zscf)
571	ENDIF
572
573	! dont let it go to low else the surface uncouples
574	q_cdrag(ji) = MAX(cd_tmp, 1.e-4/MAX(speed,min_wind))
575	!!
576	!! In some situations it might be useful to give an upper limit on the cdrag as well.
577	!! The line here should then be uncommented.
578	!! q_cdrag(ji) = MIN(q_cdrag(ji), 0.5/MAX(speed,min_wind))
579
580	END DO
581
582	IF (long_print) WRITE (numout,*) ' not ldqcdrag_from_gcm : diffuco_aero done '
583
584	END SUBROUTINE diffuco_aero
585
586	!! This routine computes partial beta coefficient : snow sublimation
587	!!
588	SUBROUTINE diffuco_snow (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, &
589	& snow, frac_nobio, totfrac_nobio, snow_nobio, vbeta1)
590
591	! interface description
592	! input scalar
593	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
594	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
595	! input fields
596	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
597	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity
598	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
599	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
600	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
601	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
602	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass
603	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice,lakes,cities,...
604	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow on ice,lakes,cities,...
605	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio!! Total fraction of ice+lakes+cities+...
606	! output fields
607	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta1 !! Beta for sublimation
608
609	! local declaration
610	REAL(r_std) :: subtest, zrapp, speed, vbeta1_add
611	INTEGER(i_std) :: ji, jv
612
613	!
614	! 1. beta coefficient for sublimation for snow on vegetation
615	!
616	DO ji=1,kjpindex
617	! Fraction of mesh that can sublimate snow
618	vbeta1(ji) = (un - totfrac_nobio(ji)) * MAX( MIN(snow(ji)/snowcri,un), zero)
619	!
620	! -- Limitation of sublimation in case of snow amounts smaller than
621	! the atmospheric demande.
622	!
623	speed = MAX(min_wind, wind(ji))
624	!
625	subtest = dtradia * vbeta1(ji) * speed * q_cdrag(ji) * rau(ji) * &
626	& ( qsatt(ji) - qair(ji) )
627	!
628	IF ( subtest .GT. zero ) THEN
629	zrapp = snow(ji) / subtest
630	IF ( zrapp .LT. un ) THEN
631	vbeta1(ji) = vbeta1(ji) * zrapp
632	ENDIF
633	ENDIF
634	!
635	END DO
636
637	!
638	! 2. add beta coefficient for other surface types.
639	!
640	DO jv = 1, nnobio
641	!!$ !
642	!!$ IF ( jv .EQ. iice ) THEN
643	!!$ !
644	!!$ ! Land ice is of course a particular case
645	!!$ !
646	!!$ DO ji=1,kjpindex
647	!!$ vbeta1(ji) = vbeta1(ji) + frac_nobio(ji,jv)
648	!!$ ENDDO
649	!!$ !
650	!!$ ELSE
651	!
652	DO ji=1,kjpindex
653	!
654	vbeta1_add = frac_nobio(ji,jv) * MAX( MIN(snow_nobio(ji,jv)/snowcri,un), zero)
655	!
656	! -- Limitation of sublimation in case of snow amounts smaller than
657	! the atmospheric demand.
658	!
659	speed = MAX(min_wind, wind(ji))
660	!
661	subtest = dtradia * vbeta1_add * speed * q_cdrag(ji) * rau(ji) * &
662	& ( qsatt(ji) - qair(ji) )
663	!
664	IF ( subtest .GT. zero ) THEN
665	zrapp = snow_nobio(ji,jv) / subtest
666	IF ( zrapp .LT. un ) THEN
667	vbeta1_add = vbeta1_add * zrapp
668	ENDIF
669	ENDIF
670	!
671	vbeta1(ji) = vbeta1(ji) + vbeta1_add
672	!
673	ENDDO
674	!!$ !
675	!!$ ENDIF
676	!
677	ENDDO
678
679	IF (long_print) WRITE (numout,*) ' diffuco_snow done '
680
681	END SUBROUTINE diffuco_snow
682
683	!! This routine computes partial beta coefficient : interception for each type of vegetation
684	!!
685	! Nathalie - Juin 2006 - Introduction de vbeta23
686	!SUBROUTINE diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
687	! & qsintveg, qsintmax, rstruct, vbeta2)
688	SUBROUTINE diffuco_inter (kjpindex, dtradia, qair, qsatt, rau, u, v, q_cdrag, veget, &
689	& qsintveg, qsintmax, rstruct, vbeta2, vbeta23)
690
691	! interface description
692	! input scalar
693	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
694	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
695	! input fields
696	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
697	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qsatt !! Surface saturated humidity
698	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
699	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
700	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
701	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
702	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! vegetation fraction for each type
703	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
704	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
705	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! STOMATE: architectural resistance
706	! output fields
707	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta2 !! Beta for interception loss
708	! AJout Nathalie - Juin 2006
709	!! Beta for fraction of wetted foliage that will transpire
710	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta23
711	! Fin ajout Nathalie
712
713	! local declaration
714	INTEGER(i_std) :: ji, jv
715	REAL(r_std) :: zqsvegrap, ziltest, zrapp, speed
716
717	!
718	! Correction Nathalie - Initialisation des vbeta2x
719	vbeta2(:,:) = zero
720	! Ajout Nathalie - Juin 2006
721	vbeta23(:,:) = zero
722	! Fin ajout Nathalie
723	!
724	DO jv = 1,nvm
725
726	!
727	! 1. beta coefficient for vegetation interception
728	!
729
730	DO ji=1,kjpindex
731
732	IF (veget(ji,jv) .GT. min_sechiba .AND. qsintveg(ji,jv) .GT. zero ) THEN
733
734	zqsvegrap = zero
735	IF (qsintmax(ji,jv) .GT. min_sechiba ) THEN
736	zqsvegrap = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
737	END IF
738	!
739	speed = MAX(min_wind, wind(ji))
740	! -- Interception loss: IL
741	vbeta2(ji,jv) = veget(ji,jv) * zqsvegrap * (un / (un + speed * q_cdrag(ji) * rstruct(ji,jv)))
742
743	!
744	! -- Limitation of IL by the water stored on the leaf.
745	! A first approximation of IL is obtained with the old values of
746	! qair and qsol_sat: function of temp-sol and pb. (see call of qsatcalc)
747	!
748	ziltest = dtradia * vbeta2(ji,jv) * speed * q_cdrag(ji) * rau(ji) * &
749	& ( qsatt(ji) - qair(ji) )
750	IF ( ziltest .GT. zero ) THEN
751	zrapp = qsintveg(ji,jv) / ziltest
752	IF ( zrapp .LT. un ) THEN
753	! Ajout Nathalie - Juin 2006
754	vbeta23(ji,jv) = MAX(vbeta2(ji,jv) - vbeta2(ji,jv) * zrapp, 0.)
755	! Fin ajout Nathalie
756	vbeta2(ji,jv) = vbeta2(ji,jv) * zrapp
757	ENDIF
758	ENDIF
759	END IF
760
761	END DO
762
763	END DO
764
765	IF (long_print) WRITE (numout,*) ' diffuco_inter done '
766
767	END SUBROUTINE diffuco_inter
768
769	!! This routine computes partial beta coefficient : bare soil
770	!!
771	SUBROUTINE diffuco_bare (kjpindex, dtradia, u, v, q_cdrag, rsol, evap_bare_lim, evapot, humrel, veget, vbeta4)
772
773	! interface description
774	! input scalar
775	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
776	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
777	! input fields
778	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
779	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
780	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
781	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rsol !! resistance for bare soil evaporation
782	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evap_bare_lim !! Beta factor for bare soil evaporation
783	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation
784	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress
785	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction
786	! output fields
787	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta4 !! Beta for bare soil evaporation
788
789	! local declaration
790	INTEGER(i_std) :: ji, jv
791	REAL(r_std) :: speed
792
793	IF ( .NOT. control%hydrol_cwrr ) THEN
794	DO ji = 1, kjpindex
795	!
796	vbeta4(ji) = zero
797	!
798	! 1. Soil resistance and beta for bare soil
799	! note: veget ( ,1) contains the fraction of bare soil
800	! see hydrol module
801	!
802	IF (veget(ji,1) .GE. min_sechiba) THEN
803	!
804	speed = MAX(min_wind, wind(ji))
805	!
806	! Correction Nathalie de Noblet - le 27 Mars 2006
807	! Selon recommandation de Frederic Hourdin: supprimer humrel dans formulation vbeta4
808	!vbeta4(ji) = veget(ji,1) humrel(ji,1) (un / (un + speed * q_cdrag(ji) * rsol(ji)))
809	! Nathalie - le 28 mars 2006 - vbeta4 n'etait pas calcule en fonction de
810	! rsol mais de rsol_cste * hdry! Dans ce cas inutile de calculer rsol(ji)!!
811	vbeta4(ji) = veget(ji,1) * (un / (un + speed * q_cdrag(ji) * rsol(ji)))
812	!
813	ENDIF
814	!
815	END DO
816	ELSE
817	DO ji = 1, kjpindex
818	vbeta4(ji) = evap_bare_lim(ji)
819	END DO
820	ENDIF
821
822	IF (long_print) WRITE (numout,*) ' diffuco_bare done '
823
824	END SUBROUTINE diffuco_bare
825
826	!! This routine computes partial beta coefficient : transpiration for each type of vegetation
827	!!
828	! Nathalie - Juin 2006 - introduction de vbeta23
829	!SUBROUTINE diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
830	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
831	SUBROUTINE diffuco_trans (kjpindex, dtradia, swnet, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
832	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
833
834	! interface description
835	! input scalar
836	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
837	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
838	! input fields
839	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Short wave net flux in
840	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin
841	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure
842	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
843	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
844	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
845	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
846	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
847	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress
848	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction
849	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. vegetation fraction (LAI -> infty)
850	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index
851	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
852	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
853	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: rstruct !! STOMATE
854	! AJout Nathalie - Juin 2006
855	!! Beta for fraction of wetted foliage that will transpire
856	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23
857	! Fin ajout Nathalie
858	! output fields
859	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration
860	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Surface resistance of vegetation
861	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE
862	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE
863
864	! local declaration
865	INTEGER(i_std) :: ji, jv
866	REAL(r_std) :: speed
867	REAL(r_std), DIMENSION(kjpindex) :: zdefconc, zqsvegrap
868	REAL(r_std), DIMENSION(kjpindex) :: qsatt
869
870	!
871	! 1. Moisture concentration at the leaf level.
872	!
873	CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
874	zdefconc(:) = rau(:) * MAX( qsatt(:) - qair(:), zero )
875
876	!
877	! 2. beta coefficient for vegetation transpiration
878	!
879
880	DO jv = 1,nvm
881
882	rveget(:,jv) = undef_sechiba
883	vbeta3(:,jv) = zero
884
885	zqsvegrap(:) = zero
886
887	DO ji = 1, kjpindex
888
889	speed = MAX(min_wind, wind(ji))
890
891	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
892	zqsvegrap(ji) = MAX(zero, qsintveg(ji,jv) / qsintmax(ji,jv))
893	ENDIF
894
895	IF ( ( veget(ji,jv)*lai(ji,jv) .GT. min_sechiba ) .AND. &
896	( kzero(jv) .GT. min_sechiba ) .AND. &
897	( swnet(ji) .GT. min_sechiba ) ) THEN
898
899	rveget(ji,jv) = (( swnet(ji) + rayt_cste ) / swnet(ji) ) &
900	* ((defc_plus + (defc_mult * zdefconc(ji) )) / kzero(jv)) * (un / lai(ji,jv))
901	! Corrections Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
902	! Introduction d'un potentiometre (rveg_pft) pour regler la somme rveg+rstruct
903	!vbeta3(ji,jv) = veget(ji,jv) * (un - zqsvegrap(ji)) * humrel(ji,jv) * &
904	! (un / (un + speed * q_cdrag(ji) * (rveget(ji,jv) + rstruct(ji,jv))))
905	vbeta3(ji,jv) = veget(ji,jv) * (un - zqsvegrap(ji)) * humrel(ji,jv) * &
906	(un / (un + speed * q_cdrag(ji) * (rveg_pft(jv)*(rveget(ji,jv) + rstruct(ji,jv)))))
907	! Fin ajout Nathalie
908	! Ajout Nathalie - Juin 2006
909	vbeta3(ji,jv) = vbeta3(ji,jv) + MIN( vbeta23(ji,jv), &
910	veget(ji,jv) * zqsvegrap(ji) * humrel(ji,jv) * &
911	(un / (un + speed * q_cdrag(ji) * (rveg_pft(jv)*(rveget(ji,jv) + rstruct(ji,jv))))))
912	! Fin ajout Nathalie
913
914	ENDIF
915
916	ENDDO
917
918	ENDDO
919
920	! STOMATE
921	cimean(:,:) = zero
922	vbetaco2(:,:) = zero
923
924	IF (long_print) WRITE (numout,*) ' diffuco_trans done '
925
926	END SUBROUTINE diffuco_trans
927
928	!! This routine computes partial beta coefficient : transpiration for each type of vegetation
929	!! STOMATE: this routine now calculates also the assimilation using the Farqhuar & al (1980) formulation
930	!!
931	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
932	! Nathalie - Juin 2006 - introduction de vbeta23
933	!SUBROUTINE diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, rau, u, v, q_cdrag, humrel, &
934	! assim_param, ccanopy, &
935	! veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2)
936	SUBROUTINE diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, &
937	assim_param, ccanopy, &
938	veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
939
940	! interface description
941	! input scalar
942	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
943	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
944	! input fields
945	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Downwelling short wave flux
946	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin
947	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure
948	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
949	! Ajout Nathalie - Juin 2006 - declaration q2m
950	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity
951	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature
952	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
953	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
954	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
955	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
956	!! min+max+opt temps, vcmax, vjmax for photosynthesis
957	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in) :: assim_param
958	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! STOMATE: CO2 concentration inside the canopy
959	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress
960	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Type of vegetation fraction
961	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. vegetation fraction (LAI -> infty)
962	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index
963	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
964	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
965	! AJout Nathalie - Juin 2006
966	!! Beta for fraction of wetted foliage that will transpire
967	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vbeta23
968	! Fin ajout Nathalie
969	! output fields
970	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbeta3 !! Beta for Transpiration
971	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rveget !! Surface resistance of vegetation
972	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: rstruct !! STOMATE
973	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: cimean !! STOMATE
974	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vbetaco2 !! STOMATE
975
976	! local declaration
977	REAL(r_std),DIMENSION (kjpindex,nvm) :: vcmax
978	REAL(r_std),DIMENSION (kjpindex,nvm) :: vjmax
979	REAL(r_std),DIMENSION (kjpindex,nvm) :: tmin
980	REAL(r_std),DIMENSION (kjpindex,nvm) :: topt
981	REAL(r_std),DIMENSION (kjpindex,nvm) :: tmax
982	INTEGER(i_std) :: ji, jv, jl
983	REAL(r_std), DIMENSION(kjpindex) :: leaf_ci_lowest
984	INTEGER(i_std), DIMENSION(kjpindex) :: ilai
985	REAL(r_std), DIMENSION(kjpindex) :: zqsvegrap
986	REAL(r_std) :: speed
987	! STOMATE:
988	LOGICAL, DIMENSION(kjpindex) :: assimilate, calculate
989	INTEGER(i_std) :: nic,inic,icinic
990	INTEGER(i_std), DIMENSION(kjpindex) :: index_calc
991	INTEGER(i_std) :: nia,inia,nina,inina,iainia
992	INTEGER(i_std), DIMENSION(kjpindex) :: index_assi,index_non_assi
993	REAL(r_std), PARAMETER :: laimax = 12.
994	REAL(r_std), PARAMETER :: xc4_1 = .83
995	REAL(r_std), PARAMETER :: xc4_2 = .93
996	REAL(r_std), DIMENSION(kjpindex) :: vc2, vj2
997	REAL(r_std), DIMENSION(kjpindex) :: assimi
998	REAL(r_std) :: x_1,x_2,x_3,x_4,x_5,x_6
999	REAL(r_std), DIMENSION(kjpindex) :: gstop, gs
1000	REAL(r_std), DIMENSION(kjpindex) :: Kc, Ko, CP
1001	REAL(r_std), DIMENSION(kjpindex) :: vc, vj
1002	REAL(r_std), DIMENSION(kjpindex) :: kt, rt
1003	REAL(r_std), DIMENSION(kjpindex) :: air_relhum
1004	REAL(r_std), DIMENSION(kjpindex) :: water_lim, temp_lim
1005	REAL(r_std), DIMENSION(kjpindex) :: gstot
1006	REAL(r_std), DIMENSION(kjpindex) :: assimtot
1007	REAL(r_std), DIMENSION(kjpindex) :: leaf_gs_top !! stomatal conductance at topmost level
1008	REAL(r_std), DIMENSION(nlai+1) :: laitab !! tabulated LAI steps
1009	REAL(r_std), DIMENSION(kjpindex) :: qsatt
1010	REAL(r_std), DIMENSION(nvm,nlai) :: light !! fraction of light that gets through
1011	REAL(r_std), DIMENSION(kjpindex) :: ci_gs
1012	REAL(r_std) :: cresist !! coefficient for resistances
1013	!
1014	! calculate LAI steps
1015	!
1016	DO jl = 1, nlai+1
1017	laitab(jl) = laimax(EXP(.15REAL(jl-1,r_std))-1.)/(EXP(.15*REAL(nlai,r_std))-1.)
1018	ENDDO
1019	!
1020	! calculate light fraction that comes through at a given LAI for each vegetation type
1021	!
1022	DO jl = 1, nlai
1023	!
1024	DO jv = 1, nvm
1025	!
1026	light(jv,jl) = exp( -ext_coef(jv)*laitab(jl) )
1027	!
1028	ENDDO
1029	!
1030	ENDDO
1031	!
1032	! 1. Photosynthesis parameters
1033	!
1034	!
1035	! temperatures in K
1036	!
1037	tmin(:,:) = assim_param(:,:,itmin)
1038	tmax(:,:) = assim_param(:,:,itmax)
1039	topt(:,:) = assim_param(:,:,itopt)
1040	!
1041	vcmax(:,:) = assim_param(:,:,ivcmax)
1042	vjmax(:,:) = assim_param(:,:,ivjmax)
1043	!
1044	! estimation of relative humidity of the air
1045	!
1046	! correction Nathalie, on utilise q2m/t2m au lieu de qair - Juin 2006
1047	! CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
1048	! air_relhum(:) = &
1049	! ( qair(:) * pb(:) / (0.622+qair(:)*0.378) ) / &
1050	! ( qsatt(:)pb(:) / (0.622+qsatt(:)0.378 ) )
1051	CALL qsatcalc (kjpindex, t2m, pb, qsatt)
1052	air_relhum(:) = &
1053	( q2m(:) * pb(:) / (0.622+q2m(:)*0.378) ) / &
1054	( qsatt(:)pb(:) / (0.622+qsatt(:)0.378 ) )
1055	!
1056	DO jv = 1,nvm
1057	!
1058	! 2. beta coefficient for vegetation transpiration
1059	!
1060	rstruct(:,jv) = rstruct_const(jv)
1061	rveget(:,jv) = undef_sechiba
1062	!
1063	vbeta3(:,jv) = zero
1064	vbetaco2(:,jv) = zero
1065	!
1066	cimean(:,jv) = ccanopy(:)
1067	!
1068	! mask that contains points where there is photosynthesis
1069	!
1070	nia=0
1071	nina=0
1072	!
1073	DO ji=1,kjpindex
1074	!
1075	IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. &
1076	( veget_max(ji,jv) .GT. 1.E-8 ) ) THEN
1077	IF ( ( veget(ji,jv) .GT. 1.E-8 ) .AND. &
1078	( swdown(ji) .GT. min_sechiba ) .AND. &
1079	( temp_air(ji) .GT. tmin(ji,jv) ) .AND. &
1080	( temp_air(ji) .LT. tmax(ji,jv) ) .AND. &
1081	( humrel(ji,jv) .GT. min_sechiba ) ) THEN
1082	!
1083	assimilate(ji) = .TRUE.
1084	nia=nia+1
1085	index_assi(nia)=ji
1086	!
1087	ELSE
1088	!
1089	assimilate(ji) = .FALSE.
1090	nina=nina+1
1091	index_non_assi(nina)=ji
1092	!
1093	ENDIF
1094	ELSE
1095	!
1096	assimilate(ji) = .FALSE.
1097	nina=nina+1
1098	index_non_assi(nina)=ji
1099	!
1100	ENDIF
1101	!
1102	ENDDO
1103	!
1104	gstot(:) = zero
1105	assimtot(:) = zero
1106	!
1107	zqsvegrap(:) = zero
1108	WHERE (qsintmax(:,jv) .GT. min_sechiba)
1109	zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
1110	ENDWHERE
1111	!
1112	WHERE ( assimilate(:) )
1113	water_lim(:) = MIN( 2.*humrel(:,jv), 1. )
1114	ENDWHERE
1115	! give a default value of ci for all pixel that do not assimilate
1116	DO jl=1,nlai
1117	DO inina=1,nina
1118	leaf_ci(index_non_assi(inina),jv,jl) = ccanopy(index_non_assi(inina)) * .667_r_std
1119	ENDDO
1120	ENDDO
1121	!
1122	ilai(:) = 1
1123	!
1124	! Here is calculated photosynthesis (Farqhuar et al. 80)
1125	! and stomatal conductance (Ball & al. 86)
1126	!
1127	! Calculating temperature dependent parameters
1128	!
1129	IF ( is_c4(jv) ) THEN
1130	!
1131	! Case of C4 plants
1132	!
1133	IF (nia .GT. 0) then
1134	!OCL NOVREC
1135	DO inia=1,nia
1136	!
1137	x_1 = 0.177 * EXP( 0.069*(temp_air(index_assi(inia))-tp_00) )
1138	! = 2.0**(((temp_air(index_assi(inia))-tp_00)-25.0)/10.0)
1139	!
1140	kt(index_assi(inia)) = 0.7 * x_1 * 1.e6
1141	rt(index_assi(inia)) = 0.8 * x_1 / &
1142	( 1.0 + EXP(1.3*(temp_air(index_assi(inia))-tmax(index_assi(inia),jv))) )
1143	!
1144	vc(index_assi(inia)) = vcmax(index_assi(inia),jv) &
1145	* 0.39 * x_1 * water_lim(index_assi(inia)) / &
1146	! * 0.39 * x_1 / &
1147	( (1.0+EXP(0.3*(tmin(index_assi(inia),jv)-temp_air(index_assi(inia))))) &
1148	* (1.0+EXP(0.3*(temp_air(index_assi(inia))-topt(index_assi(inia),jv)))) )
1149	!
1150	ENDDO
1151	ENDIF
1152	!
1153	IF (nina .GT. 0) then
1154	!
1155	!OCL NOVREC
1156	DO inina=1,nina
1157	!
1158	kt(index_non_assi(inina)) = zero
1159	rt(index_non_assi(inina)) = zero
1160	vc(index_non_assi(inina)) = zero
1161	!
1162	ENDDO
1163	!
1164	ENDIF
1165	!
1166	ELSE
1167	!
1168	! Case of C3 plants
1169	!
1170	IF (nia .GT. 0) then
1171	!
1172	!OCL NOVREC
1173	DO inia=1,nia
1174	!
1175	temp_lim(index_assi(inia)) = &
1176	(temp_air(index_assi(inia))-tmin(index_assi(inia),jv)) * &
1177	(temp_air(index_assi(inia))-tmax(index_assi(inia),jv))
1178	temp_lim(index_assi(inia)) = temp_lim(index_assi(inia)) / &
1179	(temp_lim(index_assi(inia))-(temp_air(index_assi(inia))-&
1180	topt(index_assi(inia),jv))**2)
1181	!
1182	Kc(index_assi(inia)) = 39.09 * EXP(.085*(temp_air(index_assi(inia))-tp_00))
1183	!
1184	Ko(index_assi(inia)) = 2.412 * 210000. &
1185	* EXP(.085*(temp_air(index_assi(inia))-tmin(index_assi(inia),jv))) / &
1186	(temp_air(index_assi(inia))-tmin(index_assi(inia),jv))
1187	!
1188	CP(index_assi(inia)) = 42. * EXP( 9.46*(temp_air(index_assi(inia))-tp_00-25.)/&
1189	temp_air(index_assi(inia)) )
1190	!
1191	vc(index_assi(inia)) = vcmax(index_assi(inia),jv) * &
1192	temp_lim(index_assi(inia)) * water_lim(index_assi(inia))
1193	! temp_lim(index_assi(inia))
1194	vj(index_assi(inia)) = vjmax(index_assi(inia),jv) * &
1195	temp_lim(index_assi(inia)) * water_lim(index_assi(inia))
1196	! temp_lim(index_assi(inia))
1197	!
1198	ENDDO
1199	!
1200	ENDIF
1201	!
1202	IF (nina .GT. 0) then
1203	!
1204	!OCL NOVREC
1205	DO inina=1,nina
1206	!
1207	temp_lim(index_non_assi(inina)) = zero
1208	Kc(index_non_assi(inina)) = zero
1209	Ko(index_non_assi(inina)) = zero
1210	CP(index_non_assi(inina)) = zero
1211	!
1212	vc(index_non_assi(inina)) = zero
1213	vj(index_non_assi(inina)) = zero
1214	!
1215	ENDDO
1216	!
1217	ENDIF
1218	!
1219	ENDIF ! C3/C4
1220	!
1221	! estimate assimilation and conductance for each LAI level
1222	!
1223	DO jl = 1, nlai
1224	!
1225	nic=0
1226	!
1227	calculate(:) = .FALSE.
1228	!
1229	IF (nia .GT. 0) then
1230	!
1231	!OCL NOVREC
1232	DO inia=1,nia
1233	!
1234	calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) )
1235	!
1236	IF ( calculate(index_assi(inia)) ) THEN
1237	!
1238	nic=nic+1
1239	index_calc(nic)=index_assi(inia)
1240	!
1241	ENDIF
1242	!
1243	ENDDO
1244	!
1245	ENDIF
1246	!
1247	! Vmax is scaled into the canopy due to reduction of nitrogen
1248	!
1249	x_1 = ( un - .7_r_std * ( un - light(jv,jl) ) )
1250	!
1251	IF ( nic .GT. 0 ) THEN
1252	!
1253	DO inic=1,nic
1254	!
1255	vc2(index_calc(inic)) = vc(index_calc(inic)) * x_1
1256	vj2(index_calc(inic)) = vj(index_calc(inic)) * x_1
1257	!
1258	ENDDO
1259	!
1260	ENDIF
1261	!
1262	IF ( is_c4(jv) ) THEN
1263	!
1264	! assimilation for C4 plants (Collatz & al. 91)
1265	!
1266	DO ji = 1, kjpindex
1267	!
1268	assimi(ji) = 0.
1269	!
1270	ENDDO
1271	!
1272	IF (nic .GT. 0) THEN
1273	!
1274	!OCL NOVREC
1275	DO inic=1,nic
1276	!
1277	x_1 = - ( vc2(index_calc(inic)) + 0.092 * 2.3* swdown(index_calc(inic)) * &
1278	ext_coef(jv) * light(jv,jl) )
1279	x_2 = vc2(index_calc(inic)) * 0.092 * 2.3 * swdown(index_calc(inic)) * &
1280	ext_coef(jv) * light(jv,jl)
1281	x_3 = ( -x_1 - sqrt( x_1x_1 - 4.0 xc4_1 * x_2 ) ) / (2.0*xc4_1)
1282	x_4 = - ( x_3 + kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * &
1283	1.0e-6 )
1284	x_5 = x_3 * kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * 1.0e-6
1285	assimi(index_calc(inic)) = ( -x_4 - sqrt( x_4x_4 - 4. xc4_2 * x_5 ) ) / (2.*xc4_2)
1286	assimi(index_calc(inic)) = assimi(index_calc(inic)) - &
1287	rt(index_calc(inic))
1288	!
1289	ENDDO
1290	!
1291	ENDIF
1292	!
1293	ELSE
1294	!
1295	! assimilation for C3 plants (Farqhuar & al. 80)
1296	!
1297	DO ji = 1, kjpindex
1298	!
1299	assimi(ji) = 0.
1300	!
1301	ENDDO
1302	!
1303	IF (nic .GT. 0) THEN
1304	!
1305	!OCL NOVREC
1306	DO inic=1,nic
1307	!
1308	x_1 = vc2(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) / &
1309	( leaf_ci(index_calc(inic),jv,jl) + Kc(index_calc(inic)) * &
1310	( 1._r_std + 210000._r_std / Ko(index_calc(inic)) ) )
1311	x_2 = .8855_r_stdswdown(index_calc(inic))ext_coef(jv)*light(jv,jl)
1312	x_3 = x_2+vj2(index_calc(inic))
1313	x_4 = ( x_3 - sqrt( x_3x_3 - (4._r_std.7_r_stdx_2vj2(index_calc(inic))) ) ) / &
1314	(2._r_std*.7_r_std)
1315	x_5 = x_4 * leaf_ci(index_calc(inic),jv,jl) / &
1316	( 4.5_r_std * leaf_ci(index_calc(inic),jv,jl) + &
1317	10.5_r_std*CP(index_calc(inic)) )
1318	x_6 = MIN( x_1, x_5 )
1319	assimi(index_calc(inic)) = x_6 * ( 1._r_std - CP(index_calc(inic))/&
1320	leaf_ci(index_calc(inic),jv,jl) ) - .011_r_std * vc2(index_calc(inic))
1321	!
1322	ENDDO
1323	!
1324	ENDIF
1325	!
1326	ENDIF
1327	!
1328	IF (nic .GT. 0) THEN
1329	!
1330	!OCL NOVREC
1331	!cdir NODEP
1332	DO inic=1,nic
1333	!
1334	! estimate conductance (Ball & al. 86)
1335	!
1336	icinic=index_calc(inic)
1337	! gs(icinic) = water_lim(icinic) * &
1338	gs(icinic) = &
1339	( gsslope(jv) * assimi(icinic) * &
1340	air_relhum(icinic) / ccanopy(icinic) ) &
1341	+ gsoffset(jv)
1342	gs(icinic) = MAX( gs(icinic), gsoffset(jv) )
1343	ENDDO
1344	!
1345	DO inic=1,nic
1346	icinic=index_calc(inic)
1347	!
1348	! the new ci is calculated with
1349	! dci/dt=(ccanopy-ci)*gs/1.6-A
1350	! ci=ci+((ccanopy(icinic)-ci)*gs/1.6-&
1351	! assimi(icinic))*dtradia
1352	! we verify that ci is not out of possible values
1353	!
1354	ci_gs(icinic) = MIN( ccanopy(icinic), MAX( CP(icinic), &
1355	( ccanopy(icinic) - 1.6_r_std * assimi(icinic) / &
1356	gs(icinic) ) ) ) - leaf_ci(icinic,jv,jl)
1357	ENDDO
1358	!cdir NODEP
1359	DO inic=1,nic
1360	icinic=index_calc(inic)
1361	!to avoid some problem of numerical stability, the leaf_ci is bufferized
1362	leaf_ci(icinic,jv,jl) = leaf_ci(icinic,jv,jl) + ci_gs(icinic)/6.
1363	ENDDO
1364	!
1365	DO inic=1,nic
1366	icinic=index_calc(inic)
1367	!
1368	! this might be the last level for which Ci is calculated. Store it for
1369	! initialization of the remaining levels of the Ci array.
1370	!
1371	leaf_ci_lowest(icinic) = leaf_ci(icinic,jv,jl)
1372	ENDDO
1373	!
1374	!cdir NODEP
1375	DO inic=1,nic
1376	icinic=index_calc(inic)
1377	!
1378	! total assimilation and conductance
1379	assimtot(icinic) = assimtot(icinic) + &
1380	assimi(icinic) * (laitab(jl+1)-laitab(jl))
1381	gstot(icinic) = gstot(icinic) + &
1382	gs(icinic) * (laitab(jl+1)-laitab(jl))
1383	!
1384	ilai(icinic) = jl
1385	!
1386	ENDDO
1387	!
1388	ENDIF
1389	!
1390	! keep stomatal conductance of topmost level
1391	!
1392	IF ( jl .EQ. 1 ) THEN
1393	!
1394	leaf_gs_top(:) = 0.
1395	!
1396	IF ( nic .GT. 0 ) then
1397	!
1398	!OCL NOVREC
1399	DO inic=1,nic
1400	!
1401	leaf_gs_top(index_calc(inic)) = gs(index_calc(inic))
1402	!
1403	ENDDO
1404	!
1405	ENDIF
1406	!
1407	ENDIF
1408	!
1409	IF (nia .GT. 0) THEN
1410	!
1411	!OCL NOVREC
1412	DO inia=1,nia
1413	!
1414	IF ( .NOT. calculate(index_assi(inia)) ) THEN
1415	!
1416	! a) for plants that are doing photosynthesis, but whose LAI is lower
1417	! than the present LAI step, initialize it to the Ci of the lowest
1418	! canopy level
1419	!
1420	leaf_ci(index_assi(inia),jv,jl) = leaf_ci_lowest(index_assi(inia))
1421	!
1422	ENDIF
1423	!
1424	ENDDO
1425	!
1426	ENDIF
1427	!
1428	ENDDO ! loop over LAI steps
1429	!
1430	! final calculations: resistances
1431	!
1432	IF (nia .GT. 0) THEN
1433	!
1434	!OCL NOVREC
1435	!cdir NODEP
1436	DO inia=1,nia
1437	!
1438	iainia=index_assi(inia)
1439	!
1440	! conversion from mmol/m2/s to m/s
1441	!
1442	gstot(iainia) = .0244(temp_air(iainia)/tp_00)&
1443	(1013./pb(iainia))*gstot(iainia)
1444	gstop(iainia) = .0244 * (temp_air(iainia)/tp_00)*&
1445	(1013./pb(iainia))leaf_gs_top(iainia)&
1446	laitab(ilai(iainia)+1)
1447	!
1448	rveget(iainia,jv) = 1./gstop(iainia)
1449	!
1450	ENDDO
1451	!
1452	DO inia=1,nia
1453	!
1454	iainia=index_assi(inia)
1455	!
1456	! rstruct is the difference between rtot (=1./gstot) and rveget
1457	!
1458	! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
1459	!rstruct(iainia,jv) = 1./gstot(iainia) - &
1460	! rveget(iainia,jv)
1461	rstruct(iainia,jv) = MAX( 1./gstot(iainia) - &
1462	rveget(iainia,jv), min_sechiba)
1463	!
1464	ENDDO
1465	!
1466	DO inia=1,nia
1467	!
1468	iainia=index_assi(inia)
1469	!
1470	speed = MAX(min_wind, wind(index_assi(inia)))
1471	!
1472	! beta for transpiration
1473	!
1474	! Corrections Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
1475	! Introduction d'un potentiometre (rveg_pft) pour regler la somme rveg+rstruct
1476	!vbeta3(iainia,jv) = veget_max(iainia,jv) * &
1477	! (un - zqsvegrap(iainia)) * &
1478	! (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
1479	! rstruct(iainia,jv))))
1480	cresist=(un / (un + speed * q_cdrag(iainia) * &
1481	(rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
1482
1483	vbeta3(iainia,jv) = veget_max(iainia,jv) * &
1484	(un - zqsvegrap(iainia)) * cresist + &
1485	!!$ ! Ajout Nathalie - Juin 2006
1486	!!$ vbeta3(iainia,jv) = vbeta3(iainia,jv) + &
1487	! Corrections Nathalie - le 09 novembre 2009 : veget => veget_max
1488	! MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
1489	MIN( vbeta23(iainia,jv), veget_max(iainia,jv) * &
1490	! zqsvegrap(iainia) * humrel(iainia,jv) * &
1491	zqsvegrap(iainia) * cresist )
1492	! Fin ajout Nathalie
1493	!
1494	! beta for assimilation. The difference is that surface
1495	! covered by rain (un - zqsvegrap(iainia)) is not taken into account
1496	! 1.6 is conversion for H2O to CO2 conductance
1497	! vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
1498	! (un / (un + q_cdrag(iainia) * &
1499	! (rveget(iainia,jv))))/1.6
1500	!
1501	vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
1502	(un / (un + speed * q_cdrag(iainia) * &
1503	(rveget(iainia,jv) + rstruct(iainia,jv)))) / 1.6
1504	!
1505	! cimean is the "mean ci" calculated in such a way that assimilation
1506	! calculated in enerbil is equivalent to assimtot
1507	!
1508	cimean(iainia,jv) = ccanopy(iainia) - &
1509	assimtot(iainia) / &
1510	( vbetaco2(iainia,jv)/veget_max(iainia,jv) * &
1511	rau(iainia) * speed * q_cdrag(iainia))
1512	!
1513	ENDDO
1514	!
1515	ENDIF
1516	!
1517	END DO ! loop over vegetation types
1518	!
1519	IF (long_print) WRITE (numout,*) ' diffuco_trans_co2 done '
1520
1521	END SUBROUTINE diffuco_trans_co2
1522
1523	!! This routine combines previous partial beta coeeficient and calculates
1524	!! alpha and complete beta coefficient
1525	!!
1526	! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006
1527	SUBROUTINE diffuco_comb (kjpindex, dtradia, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
1528	& snow, veget, lai, vbeta1, vbeta2, vbeta3 , vbeta4, valpha, vbeta, qsintmax)
1529
1530	! interface description
1531	! input scalar
1532	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1533	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
1534	! input fields
1535	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: rau !! Density
1536	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
1537	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
1538	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
1539	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure
1540	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
1541	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Skin temperature in Kelvin
1542	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! lower air temperature in Kelvin
1543	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass
1544	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
1545	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index
1546	! Ajout Nathalie / le 13-03-2006
1547	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation
1548	! modified fields
1549	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta1 !! Beta for sublimation
1550	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vbeta4 !! Beta for Bare soil evaporation
1551	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress
1552	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2 !! Beta for Interception for
1553	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3 !! Beta for Transpiration
1554	! output fields
1555	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: valpha !! TotalAlpha coefficient
1556	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vbeta !! Total beta coefficient
1557
1558	! local declaration
1559	INTEGER(i_std) :: ji, jv
1560	REAL(r_std) :: zevtest, zsoil_moist, zrapp
1561	REAL(r_std), DIMENSION(kjpindex) :: vbeta2sum, vbeta3sum
1562	REAL(r_std), DIMENSION(kjpindex) :: vegetsum, vegetsum2
1563	REAL(r_std), DIMENSION(kjpindex) :: qsatt
1564	LOGICAL, DIMENSION(kjpindex) :: toveg, tosnow
1565	REAL(r_std) :: coeff_dew_veg
1566
1567	vbeta2sum(:) = 0.
1568	vbeta3sum(:) = 0.
1569	DO jv = 1, nvm
1570	vbeta2sum(:) = vbeta2sum(:) + vbeta2(:,jv)
1571	vbeta3sum(:) = vbeta3sum(:) + vbeta3(:,jv)
1572	ENDDO
1573
1574	!
1575	! 1. The beta and alpha coefficients are calculated.
1576	!
1577
1578	vbeta(:) = un
1579	valpha(:) = un
1580
1581	!
1582	! 2. if snow is lower than critical value
1583	!
1584
1585	DO ji = 1, kjpindex
1586
1587	IF (snow(ji) .LT. snowcri) THEN
1588
1589	vbeta(ji) = vbeta4(ji) + vbeta2sum(ji) + vbeta3sum(ji)
1590
1591	IF (vbeta(ji) .LT. min_sechiba) THEN
1592	vbeta(ji) = zero
1593	END IF
1594
1595	END IF
1596
1597	ENDDO
1598
1599	!
1600	! 3. If we are in presence of dew.
1601	!
1602
1603	! for vectorization: some arrays
1604	vegetsum(:) = 0.
1605	DO jv = 1, nvm
1606	vegetsum(:) = vegetsum(:) + veget(:,jv)
1607	ENDDO
1608	vegetsum2(:) = 0.
1609	DO jv = 2, nvm
1610	vegetsum2(:) = vegetsum2(:) + veget(:,jv)
1611	ENDDO
1612
1613	CALL qsatcalc (kjpindex, temp_sol, pb, qsatt)
1614
1615	!
1616	! 3.1 decide where the water goes (soil, vegetation, or snow)
1617	! when air moisture exceeds saturation.
1618	!
1619	toveg(:) = .FALSE.
1620	tosnow(:) = .FALSE.
1621	DO ji = 1, kjpindex
1622	IF ( qsatt(ji) .LT. qair(ji) ) THEN
1623	IF (temp_air(ji) .GT. tp_00) THEN
1624	!
1625	! 3.1.1 If it is not freezing dew is put into the
1626	! interception reservoir and on the bare soil.
1627	toveg(ji) = .TRUE.
1628	ELSE
1629	!
1630	! 3.1.2 If it is freezing water is put into the
1631	! snow reservoir.
1632	tosnow(ji) = .TRUE.
1633	ENDIF
1634	ENDIF
1635	END DO
1636
1637	! 3.1.3 now modify valpha and vbetas where necessary.
1638	!
1639	! 3.1.3.1 Soil and snow (2d)
1640	!
1641	DO ji = 1, kjpindex
1642	IF ( toveg(ji) ) THEN
1643	vbeta1(ji) = zero
1644	vbeta4(ji) = veget(ji,1)
1645	! Correction Nathalie - le 13-03-2006: le vbeta ne sera calcule qu'une fois tous les vbeta2 redefinis
1646	!vbeta(ji) = vegetsum(ji)
1647	vbeta(ji) = vbeta4(ji)
1648	valpha(ji) = un
1649	ENDIF
1650	IF ( tosnow(ji) ) THEN
1651	vbeta1(ji) = un
1652	vbeta4(ji) = zero
1653	vbeta(ji) = un
1654	valpha(ji) = un
1655	ENDIF
1656	ENDDO
1657	!
1658	! 3.1.3.2 vegetation (3d)
1659	!
1660	DO jv = 1, nvm
1661	!
1662	DO ji = 1, kjpindex
1663	!
1664	! Correction Nathalie - 13-03-2006 / si qsintmax=0, vbeta2=0
1665	IF ( toveg(ji) ) THEN
1666	IF (qsintmax(ji,jv) .GT. min_sechiba) THEN
1667	!MM
1668	! Compute part of dew that can be intercepted by leafs.
1669	IF ( lai(ji,jv) .GT. min_sechiba) THEN
1670	IF (lai(ji,jv) .GT. 1.5) THEN
1671	coeff_dew_veg= &
1672	& 0.000017lai(ji,jv)*5 &
1673	& - 0.000824lai(ji,jv)*4 &
1674	& + 0.014843lai(ji,jv)*3 &
1675	& - 0.110112lai(ji,jv)*2 &
1676	& + 0.205673*lai(ji,jv) &
1677	& + 0.887773
1678	ELSE
1679	coeff_dew_veg=1.
1680	ENDIF
1681	ELSE
1682	coeff_dew_veg=zero
1683	ENDIF
1684	vbeta2(ji,jv) = coeff_dew_veg*veget(ji,jv)
1685	! vbeta2(ji,jv) = veget(ji,jv)
1686	ELSE
1687	vbeta2(ji,jv) = zero
1688	ENDIF
1689	vbeta(ji) = vbeta(ji) + vbeta2(ji,jv)
1690	ENDIF
1691	IF ( tosnow(ji) ) vbeta2(ji,jv) = zero
1692	!
1693	ENDDO
1694	!
1695	ENDDO
1696
1697	!
1698	! 3.2 In any case there is no transpiration when air moisture is too high.
1699	!
1700
1701	DO jv = 1, nvm
1702	DO ji = 1, kjpindex
1703	IF ( qsatt(ji) .LT. qair(ji) ) THEN
1704	vbeta3(ji,jv) = zero
1705	humrel(ji,jv) = zero
1706	ENDIF
1707	ENDDO
1708	ENDDO
1709
1710	!
1711	! 3.2_bis In any case there is no interception loss on bare soil.
1712	!
1713
1714	DO ji = 1, kjpindex
1715	IF ( qsatt(ji) .LT. qair(ji) ) THEN
1716	vbeta2(ji,1) = zero
1717	ENDIF
1718	ENDDO
1719
1720	IF (long_print) WRITE (numout,*) ' diffuco_comb done '
1721
1722	END SUBROUTINE diffuco_comb
1723
1724	SUBROUTINE diffuco_raerod (kjpindex, u, v, q_cdrag, raero)
1725	!
1726	! Simply computes the aerodynamic resistance. For the moment it is
1727	! only used as a diagnostic but that may change !
1728	!
1729	IMPLICIT NONE
1730	!
1731	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1732	! input fields
1733	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
1734	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
1735	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q_cdrag !! Surface drag
1736	! output filed
1737	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: raero !! Aerodynamic resistance
1738	!
1739	! local declaration
1740	INTEGER(i_std) :: ji
1741	REAL(r_std) :: speed
1742	!
1743	DO ji=1,kjpindex
1744	!
1745	speed = MAX(min_wind, wind(ji))
1746	raero(ji) = un / (q_cdrag(ji)*speed)
1747	!
1748	ENDDO
1749	!
1750	END SUBROUTINE diffuco_raerod
1751
1752	END MODULE diffuco

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

Download in other formats: