Context Navigation

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

Last change on this file since 105 was 105, checked in by didier.solyga, 13 years ago
Put the calls to getin RVEG and HYDROL_HUMCSTE in src_parameters
File size: 70.9 KB

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

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

Download in other formats: