Context Navigation

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

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

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

Download in other formats: