Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE/src_sechiba/diffuco.f90 @ 579

Last change on this file since 579 was 42, checked in by mmaipsl, 14 years ago

MM: Replace all 0.0 by 'zero' and 1.0 by 'un',

and all 86400. by 'one_day' in the code to reduce explicit constants.

Property svn:keywords set to Revision Date HeadURL Date Author Revision

File size: 71.1 KB

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

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

Download in other formats: