Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/diffuco.f90 @ 277

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

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

Download in other formats: