Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/diffuco.f90 @ 6268

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

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

Download in other formats: