Context Navigation

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

Last change on this file since 326 was 326, checked in by didier.solyga, 13 years ago
Synchronize intersurf.f90 and slowproc.f90 with the revisions 314 and 317 of the trunk
File size: 169.7 KB

Line
1	!
2	! Daily update of leaf area index etc.
3	!
4	!< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/slowproc.f90 $
5	!< $Date: 2011-06-21 15:28:18 +0200 (Tue, 21 Jun 2011) $
6	!< $Author: martial.mancip $
7	!< $Revision: 275 $
8	!! IPSL (2006)
9	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10	!
11	!
12	MODULE slowproc
13
14	! modules used:
15
16	USE defprec
17	USE constantes
18	USE pft_parameters
19	USE ioipsl
20	USE sechiba_io
21	USE interpol_help
22	USE stomate
23	USE stomate_data
24	USE grid
25	USE parallel
26	! USE Write_Field_p
27
28	IMPLICIT NONE
29
30	PRIVATE
31	PUBLIC slowproc_main,slowproc_clear
32
33	! To use OLD or NEW iterpollation schemes :
34	INTERFACE slowproc_interlai
35	MODULE PROCEDURE slowproc_interlai_OLD, slowproc_interlai_NEW
36	END INTERFACE
37	INTERFACE slowproc_interpol
38	MODULE PROCEDURE slowproc_interpol_OLD, slowproc_interpol_NEW
39	END INTERFACE
40	INTERFACE slowproc_interpol_g
41	MODULE PROCEDURE slowproc_interpol_OLD_g, slowproc_interpol_NEW_g
42	END INTERFACE
43
44	!
45	! variables used inside slowproc module : declaration and initialisation
46	!
47
48	LOGICAL, SAVE :: l_first_slowproc = .TRUE.!! Initialisation has to be done one time
49	REAL(r_std), SAVE :: dt_slow !! time step of slow processes and STOMATE
50	!
51	INTEGER(i_std) , SAVE :: veget_update=0 !! update frequency in years for landuse
52	!
53	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: clayfraction
54	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: laimap
55	!
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_nextyear !! next year fraction of vegetation type
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_nextyear !! next year fraction of ice+lakes+cities+...
58	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio_nextyear !! next year total fraction of ice+lakes+cities+...
59
60	CONTAINS
61
62	SUBROUTINE slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
63	ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, &
64	IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
65	t2m, t2m_min, temp_sol, stempdiag, &
66	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
67	deadleaf_cover, &
68	assim_param, &
69	lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
70	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
71	co2_flux, fco2_lu)
72
73
74	! interface description
75	! input scalar
76	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
77	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
78	INTEGER(i_std),INTENT (in) :: kjpindex !! Domain size
79	REAL(r_std),INTENT (in) :: dtradia !! Time step
80	REAL(r_std),INTENT (in) :: date0 !! Initial date
81	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
82	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
83	LOGICAL, INTENT(in) :: ldforcing_write !! Logical for _forcing_ file to write
84	LOGICAL, INTENT(in) :: ldcarbon_write !! Logical for _carbon_forcing_ file to write
85	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
86	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
87	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file file identifier
88	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file file identifier
89	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier
90	! input fields
91	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indeces of the points on the map
92	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map
93	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
94	INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land
95	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
96	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
97	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Relative humidity ("moisture stress")
98	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m !! 2 m air temperature (K)
99	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m_min !! min. 2 m air temp. during forcing time step (K)
100	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature
101	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: stempdiag !! Soil temperature
102	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: shumdiag !! Relative soil moisture
103	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: litterhumdiag !! Litter humidity
104	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
105	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
106	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: gpp !! GPP (gC/(m**2 of total ground)/time step)
107	! output scalar
108	! output fields
109	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux in gC/m**2 of average ground/second
110	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! Land Cover Change CO2 flux (gC/m**2 of average ground/s)
111	! modified scalar
112	! modified fields
113	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: lai !! Surface foliaire
114	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: height !! height (m)
115	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: veget !! Fraction of vegetation type
116	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout):: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
117	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: veget_max !! Max fraction of vegetation type
118	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
119	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(inout):: soiltype !! fraction of soil type
120	REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (inout):: assim_param!! min+max+opt temps, vcmax, vjmax for photosynthesis
121	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: deadleaf_cover !! fraction of soil covered by dead leaves
122	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout):: qsintmax !! Maximum water on vegetation for interception
123
124	! local declaration
125	INTEGER(i_std) :: j, jv
126	INTEGER(i_std), SAVE :: lcanop !! soil level used for LAI
127	REAL(r_std) :: tmp_day(1)
128	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
129
130	!
131	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_maint !! autotrophic resp. (gC/(m**2 of total ground)/time step)
132	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_hetero !! heterotrophic resp. (gC/(m**2 of total ground)/time step)
133	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_growth !! growth resp. (gC/(m**2 of total ground)/time step)
134	REAL(r_std), DIMENSION(kjpindex,nvm) :: npp !! Net Ecosystem Exchange (gC/(m**2 of total ground)/time step)
135	!
136	INTEGER(i_std) , SAVE :: veget_year !! year for landuse
137	LOGICAL :: land_use_updated
138	!
139	LOGICAL, PARAMETER :: check = .FALSE.
140
141	REAL(r_std), SAVE :: sec_old = zero
142	!
143	! do initialisation
144	!
145	IF (l_first_slowproc) THEN
146
147	!
148	! 1.1 allocation, file definitions. Restart file read for Sechiba. Set flags.
149	!
150
151	IF (long_print) WRITE (numout,*) ' l_first_slowproc : call slowproc_init '
152	CALL slowproc_init (kjit, ldrestart_read, dtradia, date0, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
153	& rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,&
154	& veget_update, veget_year)
155	!
156	! Time step in days for stomate
157	dt_days = dt_slow / one_day
158	!
159	resp_maint(:,:) = zero
160	resp_hetero(:,:) = zero
161	resp_growth(:,:) = zero
162	!
163	! 1.2 check time step
164	!
165	IF ( dt_slow .LT. dtradia ) THEN
166	WRITE(numout,*) 'slow_processes: time step smaller than forcing time step.'
167	STOP 'slowproc_main'
168	ENDIF
169	!
170	IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN
171	!
172	! 1.3 call STOMATE for initialization
173	!
174	CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, &
175	ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, &
176	IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, &
177	t2m, t2m_min, temp_sol, stempdiag, &
178	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
179	deadleaf_cover, &
180	assim_param, &
181	lai, height, veget, veget_max, qsintmax, &
182	veget_nextyear, totfrac_nobio_nextyear, &
183	hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
184	co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
185	!
186	ENDIF
187	!
188	IF ( .NOT. control%ok_stomate ) THEN
189	!
190	! 1.4 initialize some variables
191	! STOMATE diagnoses some variables for SECHIBA: height, deadleaf_cover, etc.
192	! IF SECHIBA is not coupled to STOMATE, then we must set these values otherwise.
193	!
194	CALL slowproc_derivvar (kjpindex, veget, lai, &
195	qsintmax, deadleaf_cover, assim_param, height)
196	ENDIF
197
198	RETURN
199
200	ENDIF
201	!
202	!!$ ! Land USE for next year
203	!!$ land_use_updated=.FALSE.
204	!!$ IF ( (land_use) .AND. (veget_update .GT. 0) ) THEN
205	!!$ ! if next iteration is divisibled by veget_update
206	!!$ IF ( mod(kjit+1, veget_update) .le. min_sechiba) THEN
207	!!$ !
208	!!$ veget_year=veget_year+veget_year_add
209	!!$ WRITE(numout,*) 'We are updating veget for year =' , veget_year
210	!!$ !
211	!!$ ! Save veget
212	!!$ !
213	!!$ veget_lastyear(:,:) = veget_max(:,:)
214	!!$ !
215	!!$ CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, &
216	!!$ & veget_max, frac_nobio, veget_year)
217	!!$
218	!!$ land_use_updated=.TRUE.
219	!!$ ENDIF
220	!!$ ENDIF
221	!
222	! prepares restart file for the next simulation
223	!
224	IF (ldrestart_write) THEN
225
226	IF (long_print) WRITE (numout,*) ' we have to complete restart file with SLOWPROC variables '
227
228	tmp_day(1) = day_counter
229	var_name= 'day_counter'
230	IF (is_root_prc) CALL restput (rest_id, var_name, 1 , 1 , 1, kjit, tmp_day)
231
232	var_name= 'veget'
233	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, veget, 'scatter', nbp_glo, index_g)
234	!
235	var_name= 'veget_max'
236	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, veget_max, 'scatter', nbp_glo, index_g)
237	!
238	var_name= 'lai'
239	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, lai, 'scatter', nbp_glo, index_g)
240	!
241	var_name= 'frac_nobio'
242	CALL restput_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, frac_nobio, 'scatter', nbp_glo, index_g)
243	!
244	var_name= 'soiltype_frac'
245	CALL restput_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, soiltype, 'scatter', nbp_glo, index_g)
246	!
247	var_name= 'clay_frac'
248	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, clayfraction, 'scatter', nbp_glo, index_g)
249	!
250	! The height of the vegetation could in principle be recalculated at the beginning of the run.
251	! However, this is very tedious, as many special cases have to be taken into account. This variable
252	! is therefore saved in the restart file.
253	var_name= 'height'
254	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, height, 'scatter', nbp_glo, index_g)
255	!
256	IF (read_lai) THEN
257	var_name= 'laimap'
258	CALL restput_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, laimap)
259	ENDIF
260	!
261	IF (land_use) THEN
262	tmp_day(1) = REAL(veget_year,r_std)
263	var_name='veget_year'
264	IF (is_root_prc) CALL restput (rest_id, var_name, 1 , 1 , 1, kjit, tmp_day)
265	ENDIF
266	!
267	! call STOMATE to write RESTART files
268	!
269	IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN
270	CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, &
271	ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, &
272	IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, &
273	t2m, t2m_min, temp_sol, stempdiag, &
274	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
275	deadleaf_cover, &
276	assim_param, &
277	lai, height, veget, veget_max, qsintmax, &
278	veget_nextyear, totfrac_nobio_nextyear, &
279	hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
280	co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
281	ENDIF
282
283	RETURN
284
285	END IF
286	!
287	! update day counter
288	day_counter = day_counter + dtradia
289
290	IF (check) WRITE(,) "slowproc: day_counter 3",day_counter
291	!
292	!
293	! 1. Compute date
294	!
295	! Test each day and assert all slow processes (days and years)
296	!
297	IF ( sec_old >= one_day - dtradia .AND. sec >= zero ) THEN
298	!
299	! reset counter
300	day_counter = zero
301	IF (check) WRITE(,) "slowproc: day_counter 2",day_counter
302	!
303	! Active slow processes
304	do_slow = .TRUE.
305	!
306	! count days
307	date = date + nint(dt_days)
308	IF (check) WRITE(numout,*) "New date : ",date, 'year_length ',year_length,kjit
309	!
310	! is one year over?
311	! EndOfYear must be true once per year
312	! during a call of stomate_season.
313	!
314	IF ( month == 1 .AND. day == 1 .AND. sec .LT. dtradia ) THEN
315	EndOfYear = .TRUE.
316	ELSE
317	EndOfYear = .FALSE.
318	ENDIF
319
320	ELSE
321	do_slow = .FALSE.
322	EndOfYear = .FALSE.
323	ENDIF
324	sec_old = sec
325
326	IF ( EndOfYear ) THEN
327	WRITE(numout,*) 'slowproc: EndOfYear is activated.'
328	ENDIF
329
330	! Land USE for next year
331	land_use_updated=.FALSE.
332	IF ( (land_use) .AND. (veget_update .GT. 0) ) THEN
333	! if next iteration is divisibled by veget_update
334	IF ( EndOfYear ) THEN
335	!
336	veget_year = veget_year + 1
337	!
338	IF ( MOD(veget_year - veget_year_orig, veget_update) == 0 ) THEN
339
340	WRITE(numout,*) 'We are updating land use veget for year =' , veget_year
341	!
342	CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, &
343	& veget_max, frac_nobio, veget_nextyear, frac_nobio_nextyear, veget_year)
344	!
345	! If some PFT has disapeared in new map
346	WHERE(veget_nextyear(:,:).LT.veget(:,:))
347	veget(:,:) = veget_nextyear(:,:)
348	ENDWHERE
349	!
350	DO j = 1, nnobio
351	totfrac_nobio_nextyear(:) = totfrac_nobio_nextyear(:) + frac_nobio_nextyear(:,j)
352	ENDDO
353	land_use_updated=.TRUE.
354	!
355	ENDIF
356	!
357	ENDIF
358	ENDIF ! Land Use part
359	IF (EndOfYear .AND. .NOT. land_use_updated) THEN
360	lcchange=.FALSE.
361	ENDIF
362
363	IF ( control%stomate_watchout .OR. control%ok_stomate ) THEN
364	!
365	! 2. call STOMATE, either because we want to keep track of long-term variables or
366	! because STOMATE is activated
367	!
368	CALL stomate_main (kjit, kjpij, kjpindex, dtradia, dt_slow, &
369	ldrestart_read, ldrestart_write, ldforcing_write, ldcarbon_write, &
370	IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, &
371	t2m, t2m_min, temp_sol, stempdiag, &
372	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
373	deadleaf_cover, &
374	assim_param, &
375	lai, height, veget, veget_max, qsintmax, &
376	veget_nextyear, totfrac_nobio_nextyear, &
377	hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
378	co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth)
379	IF ( control%ok_stomate .AND. control%ok_sechiba ) THEN
380	CALL histwrite(hist_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
381	CALL histwrite(hist_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg)
382	CALL histwrite(hist_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg)
383	npp(:,1)=zero
384	DO j = 2,nvm
385	npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j)
386	ENDDO
387	CALL histwrite(hist_id, 'npp', kjit, npp, kjpindex*nvm, indexveg)
388	ENDIF
389	IF ( hist2_id > 0 ) THEN
390	IF ( control%ok_stomate ) THEN
391	CALL histwrite(hist2_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
392	CALL histwrite(hist2_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg)
393	CALL histwrite(hist2_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg)
394	CALL histwrite(hist2_id, 'npp', kjit, npp, kjpindex*nvm, indexveg)
395	ENDIF
396	ENDIF
397	ENDIF
398	!
399	IF ( .NOT. control%ok_stomate ) THEN
400
401	!
402	! 2 STOMATE is not activated: we have to guess lai etc.
403	!
404
405	!
406	! 2.2 test if we have work to do
407	!
408
409	IF ( do_slow .OR. land_use_updated ) THEN
410	!
411	! 2.2.1 do daily processes if necessary
412	!
413	IF (long_print) WRITE (numout,*) 'slowproc_main : We update the daily variables'
414
415	! 2.2.2 updates lai
416	CALL slowproc_lai (kjpindex, lcanop,stempdiag, &
417	lalo,resolution,lai,month,day,read_lai,laimap)
418	!
419	IF (land_use_updated) THEN
420	veget_max(:,:)=veget_nextyear(:,:)
421	frac_nobio(:,:)=frac_nobio_nextyear(:,:)
422	ENDIF
423	!
424	! 2.2.3 updates veget
425	CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget)
426	totfrac_nobio(:) = zero
427	DO jv = 1, nnobio
428	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
429	ENDDO
430
431	! 2.2.4 updates qsintmax and other derived variables
432	CALL slowproc_derivvar (kjpindex, veget, lai, &
433	qsintmax, deadleaf_cover, assim_param, height)
434	ELSE
435	!
436	IF (land_use_updated) THEN
437	frac_nobio(:,:)=frac_nobio_nextyear(:,:)
438	ENDIF
439	!
440	END IF
441
442	!
443	! 2.3 some output fields
444	!
445
446	co2_flux(:,:) = zero
447
448	ENDIF
449
450	IF (long_print) WRITE (numout,*) ' slowproc_main done '
451
452	END SUBROUTINE slowproc_main
453	!!
454	!!
455	!!
456	SUBROUTINE slowproc_init (kjit, ldrestart_read, dtradia, date0, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
457	& rest_id, read_lai, lai, veget, frac_nobio, totfrac_nobio, soiltype, veget_max, height, lcanop,&
458	& veget_update, veget_year)
459
460	! interface description
461	! input scalar
462	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
463	LOGICAL, INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read
464	REAL(r_std),INTENT (in) :: dtradia !! Time step
465	REAL(r_std), INTENT (in) :: date0 !! intial date
466	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
467	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
468	! input fields
469	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indeces of the points on the map
470	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
471	INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land
472	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
473	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
474	! output scalar
475	INTEGER(i_std), INTENT(out) :: lcanop !! soil level used for LAI
476	INTEGER(i_std), INTENT(out) :: veget_update !! update frequency in timesteps for landuse
477	INTEGER(i_std), INTENT(out) :: veget_year !! first year for landuse
478	! output fields
479	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Surface foliere
480	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type
481	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice,lakes,cities, ...
482	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities+...
483	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Max fraction of vegetation type
484	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! Height of vegetation
485	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltype !! fraction of soil type
486	! local declaration
487	REAL(r_std) :: tmp_day(1)
488	REAL(r_std) :: tmp_veget_year(1)
489	REAL(r_std) :: zcanop !! soil depth taken for canopy
490	INTEGER(i_std) :: vtmp(1)
491	REAL(r_std), DIMENSION(nbdl) :: zsoil !! soil depths at diagnostic levels
492	INTEGER(i_std) :: j,l !! Index
493	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
494	INTEGER(i_std) :: ji, jv, ier
495	! DS 08072011 change in intent(in)
496	LOGICAL, INTENT(in) :: read_lai
497	REAL(r_std) :: frac_nobio1 !! temporary
498	REAL(r_std), DIMENSION(kjpindex,nbdl) :: stempdiag2_bid !! matrix to store stempdiag_bid
499	CHARACTER(LEN=4) :: vegsoil_dist !! Flag to choose the soil/vegetation distribution
500	!
501	CHARACTER(LEN=30), SAVE :: veget_str !! update frequency for landuse
502	REAL(r_std),DIMENSION (nbp_glo,nnobio) :: frac_nobio_g !! Fraction of ice, lakes, cities etc. in the mesh (global)
503	REAL(r_std),DIMENSION (nbp_glo,nvm) :: veget_max_g !! Fraction of vegetation type (globa)
504
505	LOGICAL, PARAMETER :: check = .FALSE.
506	!
507	! DS 15032011 add for replacing SUM(veget_max(ji,nvm-1:nvm))
508	REAL(r_std) :: sum_veget_max
509
510	!
511	! 1 allocation
512	!
513
514	ALLOCATE (clayfraction(kjpindex),stat=ier)
515	IF (ier.NE.0) THEN
516	WRITE (numout,*) ' error in clayfraction allocation. We stop. We need kjpindex words = ',kjpindex
517	STOP 'slowproc_init'
518	END IF
519	clayfraction(:)=undef_sechiba
520	!
521	veget_max(:,1) = un
522	veget_max(:,2:nvm) = zero
523	frac_nobio(:,:) = zero
524	totfrac_nobio(:) = zero
525	!
526	ier=-1
527	ALLOCATE(veget_nextyear(kjpindex, nvm), STAT=ier)
528	IF (ier/=0) THEN
529	WRITE(numout,*) "ERROR IN ALLOCATION of veget_nextyear : ",ier
530	STOP
531	ENDIF
532	veget_nextyear(:,1) = un
533	veget_nextyear(:,2:nvm) = zero
534	!
535	ier=-1
536	ALLOCATE(frac_nobio_nextyear(kjpindex, nnobio), STAT=ier)
537	IF (ier/=0) THEN
538	PRINT *,"ERROR IN ALLOCATION of frac_nobio_nextyear : ",ier
539	STOP
540	ENDIF
541	frac_nobio_nextyear(:,:) = zero
542	!
543	ier=-1
544	ALLOCATE(totfrac_nobio_nextyear(kjpindex), STAT=ier)
545	IF (ier/=0) THEN
546	PRINT *,"ERROR IN ALLOCATION of totfrac_nobio_nextyear : ",ier
547	STOP
548	ENDIF
549	!MM must be corrected when nnobio > 1 !!
550	totfrac_nobio_nextyear(:) = nnobio*un
551	!
552	! 2 read restart file
553	!
554
555	var_name= 'day_counter'
556	CALL ioconf_setatt('UNITS', 'd')
557	CALL ioconf_setatt('LONG_NAME','Fraction of computed day')
558	IF (is_root_prc) THEN
559	CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_day)
560	IF (tmp_day(1) == val_exp) THEN
561	day_counter = zero
562	ELSE
563	day_counter = tmp_day(1)
564	ENDIF
565	ENDIF
566	CALL bcast(day_counter)
567
568	! get restart value if none were found in the restart file
569	!
570	!Config Key = SECHIBA_DAY
571	!Config Desc = Time within the day simulated
572	!Config Def = 0.0
573	!Config Help = This is the time spent simulating the current day. This variable is
574	!Config prognostic as it will trigger all the computations which are
575	!Config only done once a day.
576	!
577	CALL setvar_p (day_counter, val_exp, 'SECHIBA_DAY', zero)
578	!
579	var_name= 'veget'
580	CALL ioconf_setatt('UNITS', '-')
581	CALL ioconf_setatt('LONG_NAME','Vegetation fraction')
582	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget, "gather", nbp_glo, index_g)
583	!
584	var_name= 'veget_max'
585	CALL ioconf_setatt('UNITS', '-')
586	CALL ioconf_setatt('LONG_NAME','Maximum vegetation fraction')
587	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget_max, "gather", nbp_glo, index_g)
588
589	!
590	frac_nobio(:,:) = val_exp
591	var_name= 'frac_nobio'
592	CALL ioconf_setatt('UNITS', '-')
593	CALL ioconf_setatt('LONG_NAME','Special soil type fraction')
594	CALL restget_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, .TRUE., frac_nobio, "gather", nbp_glo, index_g)
595	!
596	veget_update=0
597
598	IF (land_use) THEN
599	!
600	var_name= 'veget_year'
601	CALL ioconf_setatt('UNITS', '-')
602	CALL ioconf_setatt('LONG_NAME','Last year get in Land Use file.')
603	IF (is_root_prc) THEN
604	CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_veget_year)
605	!
606	IF (tmp_veget_year(1) == val_exp) THEN
607	veget_year=veget_year_orig
608	ELSE
609	IF (veget_reinit) THEN
610	veget_year=veget_year_orig
611	ELSE
612	veget_year=INT(tmp_veget_year(1))
613	ENDIF
614	ENDIF
615	ENDIF
616	CALL bcast(veget_year)
617	!
618	!Config Key = VEGET_UPDATE
619	!Config Desc = Update vegetation frequency
620	!Config If = LAND_USE
621	!Config Def = 0Y
622	!Config Help = The veget datas will be update each this time step.
623	!
624	veget_update=0
625	WRITE(veget_str,'(a)') '0Y'
626	CALL getin_p('VEGET_UPDATE', veget_str)
627	!
628	!
629	l=INDEX(TRIM(veget_str),'Y')
630	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
631	WRITE(numout,*) "Update frequency for land use in years :",veget_update
632	!
633	IF ( veget_update == 0 .AND. lcchange ) THEN
634	CALL ipslerr (2,'slowproc_init', &
635	& 'You have asked for LAND_COVER_CHANGE activated with VEGET_UPDATE = 0Y.',&
636	& 'We can''t use this land cover change model if veget is not updated.', &
637	& 'We have disabled it.')
638	lcchange=.FALSE.
639	ENDIF
640
641	ENDIF
642	!
643	var_name= 'soiltype_frac'
644	CALL ioconf_setatt('UNITS', '-')
645	CALL ioconf_setatt('LONG_NAME','Fraction of each soil type')
646	CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., soiltype, "gather", nbp_glo, index_g)
647	!
648	var_name= 'clay_frac'
649	CALL ioconf_setatt('UNITS', '-')
650	CALL ioconf_setatt('LONG_NAME','Fraction of clay in each mesh')
651	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., clayfraction, "gather", nbp_glo, index_g)
652	!
653	var_name= 'lai'
654	CALL ioconf_setatt('UNITS', '-')
655	CALL ioconf_setatt('LONG_NAME','Leaf area index')
656	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., lai, "gather", nbp_glo, index_g)
657	!
658	! The height of the vegetation could in principle be recalculated at the beginning of the run.
659	! However, this is very tedious, as many special cases have to be taken into account. This variable
660	! is therefore saved in the restart file.
661	var_name= 'height'
662	CALL ioconf_setatt('UNITS', 'm')
663	CALL ioconf_setatt('LONG_NAME','Height of vegetation')
664	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., height, "gather", nbp_glo, index_g)
665
666	!
667	IF (read_lai)THEN
668	!
669	ALLOCATE (laimap(kjpindex,nvm,12),stat=ier)
670	IF (ier.NE.0) THEN
671	WRITE (numout,) ' error in laimap allocation. We stop. We need kjpindexnvm12 words = ',kjpindexnvm*12
672	STOP 'slowproc_init'
673	END IF
674	laimap(:,:,:) = val_exp
675	!
676	var_name= 'laimap'
677	CALL ioconf_setatt('UNITS', '-')
678	CALL ioconf_setatt('LONG_NAME','Leaf area index read')
679	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, .TRUE., laimap)
680	!
681	ELSE
682	!
683	ALLOCATE (laimap(1,1,1))
684	ENDIF
685	!
686	!Config Key = SECHIBA_ZCANOP
687	!Config Desc = Soil level (m) used for canopy development (if STOMATE disactivated)
688	!Config Def = 0.5
689	!Config Help = The temperature at this soil depth is used to determine the LAI when
690	!Config STOMATE is not activated.
691	!
692	zcanop = 0.5_r_std
693	CALL setvar_p (zcanop, val_exp, 'SECHIBA_ZCANOP', 0.5_r_std)
694	!
695	! Initialisation of dpu
696	dpu(:)=dpu_cste
697	!MM, T. d'O. : before in constantes_soil :
698	! diaglev = &
699	! & (/ 0.001, 0.004, 0.01, 0.018, 0.045, &
700	! & 0.092, 0.187, 0.375, 0.750, 1.5, &
701	! & 2.0 /)
702	! Place here because it is used in sechiba and stomate (then in teststomate)
703	! to be in sechiba when teststomate will have disapeared.
704	!MM Problem here with dpu which depends on soil type
705	DO l = 1, nbdl-1
706	! first 2.0 is dpu
707	! second 2.0 is average
708	diaglev(l) = dpu_cste/(2*(nbdl-1) -1) ( ( 2(l-1) -1) + ( 2(l) -1) ) / 2.0
709	ENDDO
710	diaglev(nbdl) = dpu_cste
711
712	! depth at center of the levels
713	zsoil(1) = diaglev(1) / 2.
714	DO l = 2, nbdl
715	zsoil(l) = ( diaglev(l) + diaglev(l-1) ) / 2.
716	ENDDO
717
718	! index of this level
719	vtmp = MINLOC ( ABS ( zcanop - zsoil(:) ) )
720	lcanop = vtmp(1)
721
722	!
723	!!$ WRITE(numout, *)' SECHIBA_QSINT, qsintcst = ', qsintcst
724	!
725	! Time step of STOMATE and LAI update
726	!
727	!Config Key = DT_SLOW
728	!Config Desc = Time step of STOMATE and other slow processes
729	!Config Def = one_day
730	!Config Help = Time step (s) of regular update of vegetation
731	!Config cover, LAI etc. This is also the time step
732	!Config of STOMATE.
733
734	dt_slow = one_day
735	CALL getin_p('DT_SLOW', dt_slow)
736
737	!
738	! get restart value if none were found in the restart file
739	!
740	!!$ IF ( .NOT. agriculture .AND. land_use ) THEN
741	!!$ CALL ipslerr (2,'slowproc_init', &
742	!!$ & 'Problem with agriculture desactivated and Land Use activated.',&
743	!!$ & 'Are you sure ?', &
744	!!$ & '(check your parameters).')
745	!!$ ENDIF
746	!
747	IF ( impveg ) THEN
748	!
749	! We are on a point and thus we can read the information from the run.def
750	!
751	!Config Key = SECHIBA_VEG
752	!Config Desc = Vegetation distribution within the mesh (0-dim mode)
753	!Config If = IMPOSE_VEG
754	!Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0
755	!Config Help = The fraction of vegetation is read from the restart file. If
756	!Config it is not found there we will use the values provided here.
757	!
758	CALL setvar_p (veget, val_exp, 'SECHIBA_VEG', veget_ori_fixed_test_1)
759
760	!
761	!Config Key = SECHIBA_VEGMAX
762	!Config Desc = Maximum vegetation distribution within the mesh (0-dim mode)
763	!Config If = IMPOSE_VEG
764	!Config Def = 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0
765	!Config Help = The fraction of vegetation is read from the restart file. If
766	!Config it is not found there we will use the values provided here.
767	!
768	CALL setvar_p (veget_max, val_exp, 'SECHIBA_VEGMAX', veget_ori_fixed_test_1)
769
770	!
771	!Config Key = SECHIBA_FRAC_NOBIO
772	!Config Desc = Fraction of other surface types within the mesh (0-dim mode)
773	!Config If = IMPOSE_VEG
774	!Config Def = 0.0
775	!Config Help = The fraction of ice, lakes, etc. is read from the restart file. If
776	!Config it is not found there we will use the values provided here.
777	!Config For the moment, there is only ice.
778	!
779	! laisser ca tant qu'il n'y a que la glace. Pb avec setvar?
780	frac_nobio1 = frac_nobio(1,1)
781	CALL setvar_p (frac_nobio1, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1)
782	frac_nobio(:,:) = frac_nobio1
783	! CALL setvar (frac_nobio, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1)
784
785	!
786	!Config Key = SECHIBA_LAI
787	!Config Desc = LAI for all vegetation types (0-dim mode)
788	!Config Def = 0., 8., 8., 4., 4.5, 4.5, 4., 4.5, 4., 2., 2., 2., 2.
789	!Config If = IMPOSE_VEG
790	!Config Help = The maximum LAI used in the 0dim mode. The values should be found
791	!Config in the restart file. The new values of LAI will be computed anyway
792	!Config at the end of the current day. The need for this variable is caused
793	!Config by the fact that the model may stop during a day and thus we have not
794	!Config yet been through the routines which compute the new surface conditions.
795	!
796	CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax)
797
798	IF (impsoilt) THEN
799	!Config Key = SOIL_FRACTIONS
800	!Config Desc = Fraction of the 3 soil types (0-dim mode)
801	!Config Def = 0.28, 0.52, 0.20
802	!Config If = IMPOSE_VEG
803	!Config If = IMPOSE_SOILT
804	!Config Help = Determines the fraction for the 3 soil types
805	!Config in the mesh in the following order : sand loam and clay.
806	!
807	CALL setvar_p (soiltype, val_exp, 'SOIL_FRACTIONS', soiltype_default)
808
809	!Config Key = CLAY_FRACTION
810	!Config Desc = Fraction of the clay fraction (0-dim mode)
811	!Config Def = 0.2
812	!Config If = IMPOSE_VEG
813	!Config If = IMPOSE_SOILT
814	!Config Help = Determines the fraction of clay in the grid box.
815	!
816	CALL setvar_p (clayfraction, val_exp, 'CLAY_FRACTION', clayfraction_default)
817	ELSE
818	IF ( MINVAL(soiltype) .EQ. MAXVAL(soiltype) .AND. MAXVAL(soiltype) .EQ. val_exp .OR.&
819	& MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN
820
821	CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soiltype, clayfraction)
822	ENDIF
823	ENDIF
824	!
825	!Config Key = SLOWPROC_HEIGHT
826	!Config Desc = Height for all vegetation types (m)
827	!Config Def = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0
828	!Config If = IMPOSE_VEG
829	!Config Help = The height used in the 0dim mode. The values should be found
830	!Config in the restart file. The new values of height will be computed anyway
831	!Config at the end of the current day. The need for this variable is caused
832	!Config by the fact that the model may stop during a day and thus we have not
833	!Config yet been through the routines which compute the new surface conditions.
834	!
835	CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc)
836
837	ELSE
838	!
839	! We are in the full 2-D case thus we need to do the interpolation if the potential vegetation
840	! is not available
841	!
842	IF ( ALL( veget_max(:,:) .EQ. val_exp ) .OR. ALL( frac_nobio(:,:) .EQ. val_exp ) ) THEN
843
844	IF ( .NOT. land_use ) THEN
845
846	! The interpolation of vegetation has changed.
847	IF (is_root_prc) THEN
848	IF ( .NOT. old_veget ) THEN
849	! NEW slowproc interpol :
850	CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, veget_max_g, frac_nobio_g)
851	ELSE
852	! OLD slowproc interpol :
853	CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, veget_max_g, frac_nobio_g)
854	ENDIF
855	ENDIF
856	CALL scatter(veget_max_g,veget_max)
857	CALL scatter(frac_nobio_g, frac_nobio)
858	!
859	IF ( control%ok_dgvm ) THEN
860	!
861	! If we are dealing with dynamic vegetation then all
862	! natural PFTs should be set to veget_max = 0
863	! And in case no agriculture is desired, agriculture PFTS should be
864	! set to 0 as well
865
866	!$$$ 25/10/10 :: Modif DS & NV :: attention garde t'on la masse de C !!!
867	! Modif 15/03/2011
868
869	IF (agriculture) THEN
870	DO jv=2,nvm
871	IF (natural(jv)) THEN
872	veget_max(:,jv)=zero
873	ELSE
874	sum_veget_max = zero
875	DO ji = 1, kjpindex
876	! we sum only on the indexes corresponding to the non_natural pfts
877	sum_veget_max = sum_veget_max + veget_max(ji,jv)
878	! veget_max(ji,1) = un - SUM(veget_max(ji,jv))-SUM(frac_nobio(ji,:))
879	ENDDO
880	veget_max(ji,1) = un - sum_veget_max - SUM(frac_nobio(ji,:))
881	ENDIF
882	ENDDO !- end loop nvm
883	ELSE
884	veget_max(:,:) = zero
885	DO ji = 1, kjpindex
886	veget_max(ji,1) = un - SUM(frac_nobio(ji,:))
887	ENDDO
888	ENDIF
889	!
890	ENDIF
891	ELSE
892	WRITE(numout,*) 'We initialize land use veget for year =' , veget_year
893	! If restart doesn't contain veget, then it is the first computation
894	CALL slowproc_update(kjpindex, lalo, neighbours, resolution, contfrac, &
895	& veget_nextyear, frac_nobio_nextyear, veget_max, frac_nobio, veget_year, init=.TRUE.)
896	!
897	IF ( control%ok_dgvm ) THEN
898	!
899	! If we are dealing with dynamic vegetation then all
900	! natural PFTs should be set to veget_max = 0
901	! And in case no agriculture is desired, agriculture PFTS should be
902	! set to 0 as well
903
904	!$$$ 25/10/10 :: Modif DS & NV :: attention garde t'on la masse de C !!!
905	! DS : Nouvelle correction 15/03/2011
906	IF (agriculture) THEN
907	DO jv = 2, nvm
908	IF (natural(jv)) THEN
909	veget_max(:,jv)=zero
910	ELSE
911	! Add 15/03/2011
912	sum_veget_max = zero
913	DO ji = 1, kjpindex
914	sum_veget_max = sum_veget_max + veget_max(ji,jv)
915	ENDDO
916	veget_max(ji,1) = un - sum_veget_max- SUM(frac_nobio(ji,:))
917	ENDIF
918	ENDDO !- end loop nvm
919	ELSE
920	veget_max(:,:) = zero
921	DO ji = 1, kjpindex
922	veget_max(ji,1) = un - SUM(frac_nobio(ji,:))
923	ENDDO
924	ENDIF
925	!
926	ENDIF
927	!
928	ENDIF
929	!
930	ELSE
931	! WITH restarts for vegetation and DGVM and NO AGRICULTURE
932	IF ( control%ok_dgvm ) THEN
933	!
934	! If we are dealing with dynamic vegetation then all
935	! natural PFTs should be set to veget_max = 0
936	! And in case no agriculture is desired, agriculture PFTS should be
937	! set to 0 as well
938	!
939	! IF (.NOT. agriculture) THEN
940	! DO ji = 1, kjpindex
941	! veget_max(ji,1) = veget_max(ji,1) + SUM(veget_max(ji,nvm-1:nvm))
942	! ENDDO
943	! veget_max(ji,nvm-1:nvm) = zero
944	! ENDIF
945
946	!$$$ 25/10/10 :: Modif DS & NV :: attention garde t'on la masse de C !!!
947	IF (.NOT. agriculture) THEN
948	DO jv = 2, nvm
949	DO ji = 1, kjpindex
950	IF ( .NOT. natural (jv)) THEN
951	veget_max(ji,1) = veget_max(ji,1) + veget_max(ji,jv)
952	veget_max(ji,jv) = zero
953	ENDIF
954	ENDDO
955	ENDDO
956	ENDIF
957	!
958	ENDIF
959	!
960	ENDIF
961	!
962	totfrac_nobio(:) = zero
963	DO j = 1, nnobio
964	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,j)
965	ENDDO
966	!
967	!
968	IF (read_lai) THEN
969
970	!
971	! In case the LAI map was not found in the restart then we need to recompute it
972	!
973	IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN
974	! The interpolation of LAI has changed.
975	IF ( .NOT. old_lai ) THEN
976	! NEW slowproc interlai :
977	CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap)
978	ELSE
979	! OLD slowproc interlai :
980	CALL slowproc_interlai(kjpindex, lalo, resolution, laimap)
981	ENDIF
982	!
983	! Compute the current LAI just to be sure
984	!
985	stempdiag2_bid(1:kjpindex,1:nbdl) = stempdiag_bid
986	CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, &
987	lalo,resolution,lai,month,day,read_lai,laimap)
988	!
989	! Make sure that we redo the computation when we will be back in slowproc_main
990	day_counter = dt_slow - dtradia
991	!
992	ENDIF
993	!
994	ENDIF
995	!
996	IF ( MINVAL(lai) .EQ. MAXVAL(lai) .AND. MAXVAL(lai) .EQ. val_exp) THEN
997	!
998	! Get a first guess at the LAI
999	!
1000	IF ( read_lai ) THEN
1001	IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN
1002	! The interpolation of LAI has changed.
1003	IF ( .NOT. old_lai ) THEN
1004	! NEW slowproc interlai :
1005	CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap)
1006	ELSE
1007	! OLD slowproc interlai :
1008	CALL slowproc_interlai(kjpindex, lalo, resolution, laimap)
1009	ENDIF
1010	ENDIF
1011	ENDIF
1012	!
1013	stempdiag2_bid(1:kjpindex,1:nbdl) = stempdiag_bid
1014	CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, &
1015	lalo,resolution,lai,month,day,read_lai,laimap)
1016
1017	! Make sure that we redo the computation when we will be back in slowproc_main
1018	day_counter = dt_slow - dtradia
1019
1020	ENDIF
1021
1022	IF ( MINVAL(veget) .EQ. MAXVAL(veget) .AND. MAXVAL(veget) .EQ. val_exp) THEN
1023
1024	! Impose height
1025	DO jv = 1, nvm
1026	height(:,jv) = height_presc(jv)
1027	ENDDO
1028
1029	! Have a first guess at the vegetation fraction
1030	CALL slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget)
1031
1032	ENDIF
1033
1034	IF ( MINVAL(soiltype) .EQ. MAXVAL(soiltype) .AND. MAXVAL(soiltype) .EQ. val_exp .OR.&
1035	& MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN
1036
1037	CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soiltype, clayfraction)
1038
1039	ENDIF
1040
1041	ENDIF
1042	!
1043	! Select the preferences for the distribution of the 13 PFTs on the 3 soil types.
1044	!
1045	vegsoil_dist='EQUI'
1046	!
1047	SELECT CASE(vegsoil_dist)
1048	!
1049	! A reasonable choice
1050	!
1051	CASE('MAXR')
1052	pref_soil_veg(:,1) = pref_soil_veg_sand(:)
1053	pref_soil_veg(:,2) = pref_soil_veg_loan(:)
1054	pref_soil_veg(:,3) = pref_soil_veg_clay(:)
1055	!
1056	! Current default : equidistribution.
1057	!
1058	CASE('EQUI')
1059	!
1060	pref_soil_veg(:,:) = zero
1061	!
1062	CASE DEFAULT
1063	!
1064	WRITE(numout,*) 'The vegetation soil type preference you have chosen does not exist'
1065	WRITE(numout,*) 'You chose :', vegsoil_dist
1066	STOP 'slowproc_init'
1067	!
1068	END SELECT
1069	!
1070	!
1071	! calculate total fraction of the mesh that is covered by particular surface types: ice, lakes, ...
1072	!
1073	totfrac_nobio(:) = zero
1074	DO jv = 1, nnobio
1075	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
1076	ENDDO
1077
1078	l_first_slowproc = .FALSE.
1079
1080	IF (long_print) WRITE (numout,*) ' slowproc_init done '
1081
1082	END SUBROUTINE slowproc_init
1083	!!
1084	!! Clear Memory
1085	!!
1086	SUBROUTINE slowproc_clear
1087
1088
1089	l_first_slowproc = .TRUE.
1090	IF (ALLOCATED (clayfraction)) DEALLOCATE (clayfraction)
1091	IF (ALLOCATED (laimap)) DEALLOCATE (laimap)
1092	!
1093	!
1094	IF (ALLOCATED (veget_nextyear)) DEALLOCATE (veget_nextyear)
1095	IF (ALLOCATED (frac_nobio_nextyear)) DEALLOCATE (frac_nobio_nextyear)
1096	IF (ALLOCATED (totfrac_nobio_nextyear)) DEALLOCATE (totfrac_nobio_nextyear)
1097	!
1098	CALL stomate_clear
1099	!
1100	END SUBROUTINE slowproc_clear
1101	!!
1102	!! Derive some variables
1103	!!
1104	SUBROUTINE slowproc_derivvar (kjpindex, veget, lai, &
1105	qsintmax, deadleaf_cover, assim_param, height)
1106
1107	! interface description
1108	! input scalar
1109	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
1110	! input fields
1111	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
1112	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Surface foliere
1113	! output scalar
1114	! output fields
1115	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water on vegetation for interception
1116	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover!! fraction of soil covered by dead leaves
1117	REAL(r_std), DIMENSION (kjpindex,nvm,npco2), INTENT (out) :: assim_param !! min+max+opt temps & vmax for photosynthesis
1118	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height (m)
1119	!
1120	! local declaration
1121	!
1122	INTEGER(i_std) :: ji, jv
1123
1124	!
1125	! 1 Assimilation parameters
1126	!
1127	DO jv = 1, nvm
1128	assim_param(:,jv,itmin) = co2_tmin_fix(jv) + tp_00
1129	assim_param(:,jv,itopt) = co2_topt_fix(jv) + tp_00
1130	assim_param(:,jv,itmax) = co2_tmax_fix(jv) + tp_00
1131	assim_param(:,jv,ivcmax) = vcmax_fix(jv)
1132	assim_param(:,jv,ivjmax) = vjmax_fix(jv)
1133	ENDDO
1134
1135	!
1136	! 2 fraction of soil covered by dead leaves
1137	!
1138	deadleaf_cover(:) = zero
1139
1140	!
1141	! 3 height
1142	!
1143	DO jv = 1, nvm
1144	height(:,jv) = height_presc(jv)
1145	ENDDO
1146	!
1147	! 4 qsintmax
1148	!
1149	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
1150	! Ajout Nathalie - Juillet 2006
1151	qsintmax(:,1) = zero
1152
1153	END SUBROUTINE slowproc_derivvar
1154	!!
1155	!!
1156	!!
1157	SUBROUTINE slowproc_mean (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_mean)
1158
1159	! Accumulates field_in over a period of dt_tot.
1160	! Has to be called at every time step (dt).
1161	! Mean value is calculated if ldmean=.TRUE.
1162	! field_mean must be initialized outside of this routine!
1163
1164	!
1165	! 0 declarations
1166	!
1167
1168	! 0.1 input
1169
1170	! Domain size
1171	INTEGER(i_std), INTENT(in) :: npts
1172	! 2nd dimension (1 or npft)
1173	INTEGER(i_std), INTENT(in) :: n_dim2
1174	! Time step of STOMATE (days)
1175	REAL(r_std), INTENT(in) :: dt_tot
1176	! Time step in days
1177	REAL(r_std), INTENT(in) :: dt
1178	! Calculate mean ?
1179	LOGICAL, INTENT(in) :: ldmean
1180	! Daily field
1181	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_in
1182
1183	! 0.2 modified field
1184
1185	! Annual field
1186	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_mean
1187
1188	! =========================================================================
1189
1190	!
1191	! 1 accumulation
1192	!
1193
1194	field_mean(:,:) = field_mean(:,:) + field_in(:,:) * dt
1195
1196	!
1197	! 2 mean fields
1198	!
1199
1200	IF (ldmean) THEN
1201	field_mean(:,:) = field_mean(:,:) / dt_tot
1202	ENDIF
1203
1204	END SUBROUTINE slowproc_mean
1205	!!
1206	!!
1207	!!
1208	SUBROUTINE slowproc_long (npts, n_dim2, dt, tau, field_inst, field_long)
1209
1210	! Calculates a temporally smoothed field (field_long) from instantaneous
1211	! input fields.
1212	! Time constant tau determines the strength of the smoothing.
1213	! For tau -> infty, field_long becomes the true mean value of field_inst (but
1214	! the spinup becomes infinietly long, too).
1215	! field_long must be initialized outside of this routine!
1216
1217	!
1218	! 0 declarations
1219	!
1220
1221	! 0.1 input
1222
1223	! Domain size
1224	INTEGER(i_std), INTENT(in) :: npts
1225	! 2nd dimension (1 or npft)
1226	INTEGER(i_std), INTENT(in) :: n_dim2
1227	! Time step
1228	REAL(r_std), INTENT(in) :: dt
1229	! Integration time constant (has to have same unit as dt!)
1230	REAL(r_std), INTENT(in) :: tau
1231	! Instantaneous field
1232	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_inst
1233
1234	! 0.2 modified field
1235
1236	! Long-term field
1237	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_long
1238
1239	! =========================================================================
1240
1241	!
1242	! 1 test coherence
1243	!
1244
1245	IF ( ( tau .LT. dt ) .OR. ( dt .LE. zero ) .OR. ( tau .LE. zero ) ) THEN
1246	WRITE(numout,*) 'slowproc_long: Problem with time steps'
1247	WRITE(numout,*) 'dt=',dt
1248	WRITE(numout,*) 'tau=',tau
1249	ENDIF
1250
1251	!
1252	! 2 integration
1253	!
1254
1255	field_long(:,:) = ( field_inst(:,:)dt + field_long(:,:)(tau-dt) ) / tau
1256
1257	END SUBROUTINE slowproc_long
1258	!!
1259	!!
1260	!!
1261	SUBROUTINE slowproc_veget (kjpindex, lai, frac_nobio, veget_max, veget)
1262	!
1263	! 0. Declarations
1264	!
1265
1266	! 0.1 Input
1267	INTEGER(i_std), INTENT(in) :: kjpindex
1268	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai
1269
1270	! 0.2 Modified
1271	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio
1272	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max
1273
1274	! 0.3 Output
1275	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget
1276
1277	! 0.4 Local
1278	REAL(r_std), DIMENSION(kjpindex) :: fracsum
1279	INTEGER(i_std) :: nbad
1280	INTEGER(i_std) :: ji, jv
1281	! Test Nathalie
1282	REAL(r_std) :: SUMveg
1283	!!$ REAL(r_std) :: trans_veg
1284
1285	!
1286	! 1. Sum up veget_max and frac_nobio and test if sum is equal to 1
1287	!
1288	!
1289	! 1.1 Sum up
1290	!
1291	fracsum(:) = zero
1292	DO jv = 1, nnobio
1293	DO ji = 1, kjpindex
1294	fracsum(ji) = fracsum(ji) + frac_nobio(ji,jv)
1295	ENDDO
1296	ENDDO
1297	DO jv = 1, nvm
1298	DO ji = 1, kjpindex
1299	fracsum(ji) = fracsum(ji) + veget_max(ji,jv)
1300	ENDDO
1301	ENDDO
1302	!
1303	! 1.2 stop if there is a severe problem, that is an error above 0.01%
1304	! The rest will be scaled
1305	!
1306	nbad = COUNT( ABS(fracsum(:)-un) .GT. 0.0001 )
1307	IF ( nbad .GT. 0 ) THEN
1308	WRITE(numout,*) 'Problem with total surface areas.'
1309	DO ji = 1, kjpindex
1310	IF ( ABS(fracsum(ji)-un) .GT. 0.0001 ) THEN
1311	WRITE(numout,*) 'Point :', ji
1312	WRITE(numout,*) ' frac_nobio :', frac_nobio(ji,:)
1313	WRITE(numout,*) ' Veget_max :', veget_max(ji,:)
1314	WRITE(numout,*) ' Fracsum :', fracsum(ji), EPSILON(un)
1315	WRITE(numout,*) ' The error is :', un - ( SUM(frac_nobio(ji,:)) + SUM(veget_max(ji,:)) )
1316	STOP 'slowproc_veget'
1317	ENDIF
1318	ENDDO
1319	ENDIF
1320	!
1321	! 1.3 correction at places where the problem is surely precision-related
1322	!
1323	nbad = COUNT( ABS(fracsum(:)-un) .GT. EPSILON(un) )
1324	!
1325	IF ( nbad .GT. 0 ) THEN
1326	!
1327	IF ( long_print ) THEN
1328	WRITE(numout,*) 'WARNING! scaling frac_nobio and veget_max at', nbad, ' points'
1329	ENDIF
1330	!
1331	DO jv = 1, nnobio
1332	DO ji = 1, kjpindex
1333	IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN
1334	frac_nobio(ji,jv) = frac_nobio(ji,jv)/fracsum(ji)
1335	ENDIF
1336	ENDDO
1337	ENDDO
1338	!
1339	DO jv = 1, nvm
1340	DO ji = 1, kjpindex
1341	IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN
1342	veget_max(ji,jv) = veget_max(ji,jv)/fracsum(ji)
1343	ENDIF
1344	ENDDO
1345	ENDDO
1346	!
1347	ENDIF
1348
1349	!
1350	! 2. Set veget equal to veget_max
1351	!
1352	DO jv = 1, nvm
1353	DO ji = 1, kjpindex
1354	veget(ji,jv) = veget_max(ji,jv)
1355	ENDDO
1356	ENDDO
1357	!
1358	! 3. if lai of a vegetation type (jv > 1) is small, increase soil part
1359	!
1360	! Nathalie - nouveau calcul de veget
1361	!!$ DO jv = 1,nvm
1362	!!$ DO ji = 1, kjpindex
1363	!!$
1364	!!$ IF ((jv .GT. 1) .AND. (lai(ji,jv) .LT. 0.5)) THEN
1365	!!$ trans_veg = 2.(0.5 - lai(ji,jv)) veget(ji,jv)
1366	!!$ ! We limit the smallest fraction to 0.5%
1367	!!$ IF ( veget(ji,jv) - trans_veg .GT. min_vegfrac ) THEN
1368	!!$ veget(ji,1) = veget(ji,1) + trans_veg
1369	!!$ veget(ji,jv) = veget(ji,jv) - trans_veg
1370	!!$ ELSE
1371	!!$ veget(ji,1) = veget(ji,1) + veget(ji,jv)
1372	!!$ veget(ji,jv) = zero
1373	!!$ ENDIF
1374	!!$ ENDIF
1375	!!$
1376	!!$ ENDDO
1377	!!$ ENDDO
1378	! Ajout Nouveau calcul (stomate-like)
1379	DO ji = 1, kjpindex
1380	SUMveg = zero
1381	veget(ji,1) = veget_max(ji,1)
1382	DO jv = 2, nvm
1383	veget(ji,jv) = veget_max(ji,jv) * ( un - exp( - lai(ji,jv) * ext_coeff(jv) ) )
1384	veget(ji,1) = veget(ji,1) + (veget_max(ji,jv) - veget(ji,jv))
1385	SUMveg = SUMveg + veget(ji,jv)
1386	ENDDO
1387	SUMveg = SUMveg + veget(ji,1) + SUM(frac_nobio(ji,:))
1388	IF (SUMveg .LT. 0.99999) THEN
1389	WRITE(numout,*)' ATTENTION, en ji, SUMveg LT 1: ', ji, SUMveg
1390	WRITE(numout,*)' frac_nobio = ',SUM(frac_nobio(ji,:))
1391	WRITE(numout,*) ' ',veget(ji,:)
1392	ENDIF
1393	ENDDO
1394	!
1395	! 4. Sum up surface fractions and test if the sum is equal to 1
1396	!
1397
1398	!
1399	! 4.1 Sum up
1400	!
1401	fracsum(:) = zero
1402	DO jv = 1, nnobio
1403	DO ji = 1, kjpindex
1404	fracsum(ji) = fracsum(ji) + frac_nobio(ji,jv)
1405	ENDDO
1406	ENDDO
1407	DO jv = 1, nvm
1408	DO ji = 1, kjpindex
1409	fracsum(ji) = fracsum(ji) + veget_max(ji,jv)
1410	ENDDO
1411	ENDDO
1412
1413	!
1414	! 4.2 stop if there is a severe problem
1415	!
1416	nbad = COUNT( ABS(fracsum(:)-un) .GT. (REAL(nvm+nnobio,r_std)*EPSILON(un)) )
1417	!
1418	IF ( nbad .GT. 0 ) THEN
1419	WRITE(numout,*) 'Problem with veget or frac_nobio.'
1420	DO ji = 1, kjpindex
1421	IF ( ABS(fracsum(ji)-un) .GT. (10.*EPSILON(un)) ) THEN
1422	WRITE(numout,*) 'Point :', ji
1423	WRITE(numout,*) ' frac_nobio :', frac_nobio(ji,:)
1424	WRITE(numout,*) ' Veget :', veget(ji,:)
1425	WRITE(numout,*) ' The error is :', un - (SUM(frac_nobio(ji,:)) + SUM(veget(ji,:)))
1426	STOP 'slowproc_veget'
1427	ENDIF
1428	ENDDO
1429	ENDIF
1430
1431	!
1432	! 4.3 correction at places where the problem is surely precision-related
1433	!
1434	nbad = COUNT( ABS(fracsum(:)-un) .GT. EPSILON(un) )
1435	!
1436	IF ( nbad .GT. 0 ) THEN
1437	!
1438	IF ( long_print ) THEN
1439	WRITE(numout,*) 'slowproc_veget: scaling veget at', nbad, ' points'
1440	ENDIF
1441	!
1442	DO jv = 1, nvm
1443	DO ji = 1, kjpindex
1444	IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN
1445	veget(ji,jv) = veget(ji,jv) / fracsum(ji)
1446	ENDIF
1447	ENDDO
1448	ENDDO
1449	!
1450	DO jv = 1, nnobio
1451	DO ji = 1, kjpindex
1452	IF ( ABS(fracsum(ji)-un) .GT. EPSILON(un) ) THEN
1453	frac_nobio(ji,jv) = frac_nobio(ji,jv) / fracsum(ji)
1454	ENDIF
1455	ENDDO
1456	ENDDO
1457	!
1458	ENDIF
1459
1460	END SUBROUTINE slowproc_veget
1461	!!
1462	!!
1463	!!
1464	SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,mm,dd,read_lai,laimap)
1465	!
1466	! 0. Declarations
1467	!
1468
1469	! 0.1 Input
1470	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1471	INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI
1472	INTEGER(i_std), INTENT(in) :: mm, dd
1473	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (in) :: stempdiag !! Soil temperature
1474	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
1475	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
1476
1477	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: laimap !! LAI lue
1478	LOGICAL, INTENT(in) :: read_lai
1479	! 0.2 Output
1480	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! LAI
1481
1482	! 0.3 Local
1483	INTEGER(i_std) :: ji,jv
1484
1485
1486	! Test Nathalie. On impose LAI PFT 1 a 0
1487	! On boucle sur 2,nvm au lieu de 1,nvm
1488	lai(: ,1) = zero
1489	DO jv = 2,nvm
1490	!!$ DO jv = 1,nvm
1491
1492	SELECT CASE (type_of_lai(jv))
1493
1494	CASE ("mean ")
1495	!
1496	! 1. do the interpolation between laimax and laimin
1497	!
1498	IF ( .NOT. read_lai ) THEN
1499	lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv))
1500	ELSE
1501	DO ji = 1,kjpindex
1502	lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
1503	ENDDO
1504	ENDIF
1505	!
1506	CASE ("inter")
1507	!
1508	! 2. do the interpolation between laimax and laimin
1509	!
1510	IF ( .NOT. read_lai ) THEN
1511	DO ji = 1,kjpindex
1512	lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv))
1513	ENDDO
1514	ELSE
1515	IF (mm .EQ. 1 ) THEN
1516	IF (dd .LE. 15) THEN
1517	lai(:,jv) = laimap(:,jv,12)(1-(dd+15)/30.) + laimap(:,jv,1)((dd+15)/30.)
1518	ELSE
1519	lai(:,jv) = laimap(:,jv,1)(1-(dd-15)/30.) + laimap(:,jv,2)((dd-15)/30.)
1520	ENDIF
1521	ELSE IF (mm .EQ. 12) THEN
1522	IF (dd .LE. 15) THEN
1523	lai(:,jv) = laimap(:,jv,11)(1-(dd+15)/30.) + laimap(:,jv,12)((dd+15)/30.)
1524	ELSE
1525	lai(:,jv) = laimap(:,jv,12)(1-(dd-15)/30.) + laimap(:,jv,1)((dd-15)/30.)
1526	ENDIF
1527	ELSE
1528	IF (dd .LE. 15) THEN
1529	lai(:,jv) = laimap(:,jv,mm-1)(1-(dd+15)/30.) + laimap(:,jv,mm)((dd+15)/30.)
1530	ELSE
1531	lai(:,jv) = laimap(:,jv,mm)(1-(dd-15)/30.) + laimap(:,jv,mm+1)((dd-15)/30.)
1532	ENDIF
1533	ENDIF
1534	ENDIF
1535	!
1536	CASE default
1537	!
1538	! 3. Problem
1539	!
1540	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
1541	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
1542	STOP 'slowproc_lai'
1543
1544	END SELECT
1545
1546	ENDDO
1547
1548	END SUBROUTINE slowproc_lai
1549	!!
1550	!! Interpolate the LAI map to the grid of the model
1551	!MM TAG 1.6 model !
1552	!!
1553	SUBROUTINE slowproc_interlai_OLD(nbpt, lalo, resolution, laimap)
1554	!
1555	!
1556	!
1557	! 0.1 INPUT
1558	!
1559	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
1560	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
1561	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
1562	!
1563	! 0.2 OUTPUT
1564	!
1565	REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned
1566	!
1567	! 0.3 LOCAL
1568	!
1569	!
1570	CHARACTER(LEN=80) :: filename
1571	INTEGER(i_std) :: iml, jml, ijml, i, j, ik, lml, tml, fid, ib, jb,ip, jp, vid, ai,iki,jkj
1572	REAL(r_std) :: lev(1), date, dt, coslat
1573	INTEGER(i_std) :: itau(1)
1574	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu
1575	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu, mask
1576	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful
1577	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: laimaporig
1578	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu
1579	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
1580	INTEGER, DIMENSION(nbpt) :: n_origlai
1581	INTEGER, DIMENSION(nbpt) :: n_found
1582	REAL(r_std), DIMENSION(nbpt,nvm) :: frac_origlai
1583
1584	CHARACTER(LEN=80) :: meter
1585	REAL(r_std) :: prog, sumf
1586	LOGICAL :: found
1587	INTEGER :: idi,jdi, ilast, jlast, jj, ii, jv, inear, iprog
1588	REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
1589	!
1590	!
1591	!Config Key = LAI_FILE
1592	!Config Desc = Name of file from which the vegetation map is to be read
1593	!Config If = !LAI_MAP
1594	!Config Def = ../surfmap/lai2D.nc
1595	!Config Help = The name of the file to be opened to read the LAI
1596	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
1597	!Config map which is derived from a Nicolas VIOVY one.
1598	!
1599	filename = 'lai2D.nc'
1600	CALL getin_p('LAI_FILE',filename)
1601	!
1602	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
1603	CALL bcast(iml)
1604	CALL bcast(jml)
1605	CALL bcast(lml)
1606	CALL bcast(tml)
1607	!
1608	!
1609	ALLOCATE(lon_lu(iml))
1610	ALLOCATE(lat_lu(jml))
1611	ALLOCATE(laimap_lu(iml,jml,nvm,tml))
1612	ALLOCATE(mask_lu(iml,jml))
1613	!
1614	WRITE(numout,*) 'slowproc_interlai : Reading the LAI file'
1615	!
1616	IF (is_root_prc) THEN
1617	CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu)
1618	CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu)
1619	CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu)
1620	CALL flinget(fid, 'mask', iml, jml, 0, 0, 0, 1, mask_lu)
1621	!
1622	CALL flinclo(fid)
1623	ENDIF
1624	CALL bcast(lon_lu)
1625	CALL bcast(lat_lu)
1626	CALL bcast(laimap_lu)
1627	CALL bcast(mask_lu)
1628	!
1629	WRITE(numout,*) 'slowproc_interlai : ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml)
1630	!
1631	!
1632	ijml=iml*jml
1633	ALLOCATE(lon_ful(ijml))
1634	ALLOCATE(lat_ful(ijml))
1635	ALLOCATE(laimaporig(ijml,nvm,tml))
1636	ALLOCATE(mask(ijml))
1637
1638	DO i=1,iml
1639	DO j=1,jml
1640	iki=(j-1)*iml+i
1641	lon_ful(iki)=lon_lu(i)
1642	lat_ful(iki)=lat_lu(j)
1643	laimaporig(iki,:,:)=laimap_lu(i,j,:,:)
1644	mask(iki)=mask_lu(i,j)
1645	ENDDO
1646	ENDDO
1647	!
1648	WHERE ( laimaporig(:,:,:) .LT. 0 )
1649	laimaporig(:,:,:) = zero
1650	ENDWHERE
1651	!
1652	!
1653	ALLOCATE(lon_up(nbpt))
1654	ALLOCATE(lon_low(nbpt))
1655	ALLOCATE(lat_up(nbpt))
1656	ALLOCATE(lat_low(nbpt))
1657	!
1658	DO ib =1, nbpt
1659	!
1660	! We find the 4 limits of the grid-box. As we transform the resolution of the model
1661	! into longitudes and latitudes we do not have the problem of periodicity.
1662	! coslat is a help variable here !
1663	!
1664	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
1665	!
1666	lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
1667	lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
1668	!
1669	coslat = pi/180. * R_Earth
1670	!
1671	lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
1672	lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
1673	!
1674	!
1675	!
1676	ENDDO
1677	lon_up = NINT( lon_up * 1E6 ) * 1E-6
1678	lon_low = NINT( lon_low * 1E6 ) * 1E-6
1679	lat_up = NINT( lat_up * 1E6 ) * 1E-6
1680	lat_low = NINT( lat_low * 1E6 ) * 1E-6
1681	!
1682	! Get the limits of the integration domaine so that we can speed up the calculations
1683	!
1684	domaine_lon_min = MINVAL(lon_low)
1685	domaine_lon_max = MAXVAL(lon_up)
1686	domaine_lat_min = MINVAL(lat_low)
1687	domaine_lat_max = MAXVAL(lat_up)
1688	!
1689	!!$ WRITE(,) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max
1690	!!$ WRITE(,) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max
1691	!
1692	! Ensure that the fine grid covers the whole domain
1693	WHERE ( lon_ful(:) .LT. domaine_lon_min )
1694	lon_ful(:) = lon_ful(:) + 360.
1695	ENDWHERE
1696	!
1697	WHERE ( lon_ful(:) .GT. domaine_lon_max )
1698	lon_ful(:) = lon_ful(:) - 360.
1699	ENDWHERE
1700	lon_ful = NINT( lon_ful * 1E6 ) * 1E-6
1701	lat_ful = NINT( lat_ful * 1E6 ) * 1E-6
1702	!
1703	WRITE(numout,*) 'Interpolating the lai map :'
1704	WRITE(numout,'(2a40)')'0%--------------------------------------', &
1705	& '------------------------------------100%'
1706	!
1707	ilast = 1
1708	n_origlai(:) = 0
1709	laimap(:,:,:) = zero
1710	!
1711	DO ip=1,ijml
1712	!
1713	! Give a progress meter
1714	!
1715	! prog = ip/float(ijml)*79.
1716	! IF ( ABS(prog - NINT(prog)) .LT. 1/float(ijml)*79. ) THEN
1717	! meter(NINT(prog)+1:NINT(prog)+1) = 'x'
1718	! WRITE(numout, advance="no", FMT='(a)') ACHAR(13)
1719	! WRITE(numout, advance="no", FMT='(a80)') meter
1720	! ENDIF
1721	iprog = NINT(float(ip)/float(ijml)79.) - NINT(float(ip-1)/float(ijml)79.)
1722	IF ( iprog .NE. 0 ) THEN
1723	WRITE(numout,'(a1,$)') 'y'
1724	ENDIF
1725	!
1726	! Only start looking for its place in the smaler grid if we are within the domaine
1727	! That should speed up things !
1728	!
1729	IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
1730	( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
1731	( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
1732	( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
1733	!
1734	! look for point on GCM grid which this point on fine grid belongs to.
1735	! First look at the point on the model grid where we arrived just before. There is
1736	! a good chace that neighbouring points on the fine grid fall into the same model
1737	! grid box.
1738	!
1739	IF ( ( lon_ful(ip) .GE. lon_low(ilast) ) .AND. &
1740	( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
1741	( lat_ful(ip) .GE. lat_low(ilast) ) .AND. &
1742	( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
1743	!
1744	! We were lucky
1745	!
1746	IF (mask(ip) .GT. 0) THEN
1747	n_origlai(ilast) = n_origlai(ilast) + 1
1748	DO i=1,nvm
1749	DO j=1,12
1750	laimap(ilast,i,j) = laimap(ilast,i,j) + laimaporig(ip,i,j)
1751	ENDDO
1752	ENDDO
1753	ENDIF
1754	!
1755	ELSE
1756	!
1757	! Otherwise, look everywhere.
1758	! Begin close to last grid point.
1759	!
1760	found = .FALSE.
1761	idi = 1
1762	!
1763	DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
1764
1765	!
1766	! forward and backward
1767	!
1768	DO ii = 1,2
1769	!
1770	IF ( ii .EQ. 1 ) THEN
1771	ib = ilast - idi
1772	ELSE
1773	ib = ilast + idi
1774	ENDIF
1775	!
1776	IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
1777	IF ( ( lon_ful(ip) .GE. lon_low(ib) ) .AND. &
1778	( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
1779	( lat_ful(ip) .GE. lat_low(ib) ) .AND. &
1780	( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
1781	!
1782	IF (mask(ip) .gt. 0) THEN
1783	DO i=1,nvm
1784	DO j=1,12
1785	laimap(ib,i,j) = laimap(ib,i,j) + laimaporig(ip,i,j)
1786	ENDDO
1787	ENDDO
1788	n_origlai(ib) = n_origlai(ib) + 1
1789	ENDIF
1790	ilast = ib
1791	found = .TRUE.
1792	!
1793	ENDIF
1794	ENDIF
1795	!
1796	ENDDO
1797	!
1798	idi = idi + 1
1799	!
1800	ENDDO
1801	!
1802	ENDIF ! lucky/not lucky
1803	!
1804	ENDIF ! in the domain
1805	ENDDO
1806
1807
1808	! determine fraction of LAI points in each box of the coarse grid
1809	DO ip=1,nbpt
1810	IF ( n_origlai(ip) .GT. 0 ) THEN
1811	DO jv =1,nvm
1812	laimap(ip,jv,:) = laimap(ip,jv,:)/REAL(n_origlai(ip),r_std)
1813	ENDDO
1814	ELSE
1815	!
1816	! Now we need to handle some exceptions
1817	!
1818	IF ( lalo(ip,1) .LT. -56.0) THEN
1819	! Antartica
1820	DO jv =1,nvm
1821	laimap(ip,jv,:) = zero
1822	ENDDO
1823	!
1824	ELSE IF ( lalo(ip,1) .GT. 70.0) THEN
1825	! Artica
1826	DO jv =1,nvm
1827	laimap(ip,jv,:) = zero
1828	ENDDO
1829	!
1830	ELSE IF ( lalo(ip,1) .GT. 55.0 .AND. lalo(ip,2) .GT. -65.0 .AND. lalo(ip,2) .LT. -20.0) THEN
1831	! Greenland
1832	DO jv =1,nvm
1833	laimap(ip,jv,:) = zero
1834	ENDDO
1835	!
1836	ELSE
1837	!
1838	WRITE(numout,*) 'PROBLEM, no point in the lai map found for this grid box'
1839	WRITE(numout,*) 'Longitude range : ', lon_low(ip), lon_up(ip)
1840	WRITE(numout,*) 'Latitude range : ', lat_low(ip), lat_up(ip)
1841	!
1842	WRITE(numout,*) 'Looking for nearest point on the lai map file'
1843	CALL slowproc_nearest (ijml, lon_ful, lat_ful, &
1844	lalo(ip,2), lalo(ip,1), inear)
1845	WRITE(numout,*) 'Coordinates of the nearest point, ',inear,' :', &
1846	lon_ful(inear),lat_ful(inear)
1847	!
1848	DO jv =1,nvm
1849	laimap(ip,jv,:) = laimaporig(inear,jv,:)
1850	ENDDO
1851	ENDIF
1852	ENDIF
1853	ENDDO
1854	!
1855	WRITE(numout,*) 'slowproc_interlai : Interpolation Done'
1856	!
1857	!
1858	!
1859	DEALLOCATE(lon_up)
1860	DEALLOCATE(lon_low)
1861	DEALLOCATE(lat_up)
1862	DEALLOCATE(lat_low)
1863	DEALLOCATE(lat_ful)
1864	DEALLOCATE(lon_ful)
1865	DEALLOCATE(lat_lu)
1866	DEALLOCATE(lon_lu)
1867	DEALLOCATE(laimap_lu)
1868	DEALLOCATE(laimaporig)
1869	DEALLOCATE(mask_lu)
1870	DEALLOCATE(mask)
1871	!
1872	RETURN
1873	!
1874	END SUBROUTINE slowproc_interlai_OLD
1875	!!
1876	!! Interpolate the LAI map to the grid of the model
1877	!!
1878	SUBROUTINE slowproc_interlai_NEW(nbpt, lalo, resolution, neighbours, contfrac, laimap)
1879	!
1880	!
1881	!
1882	! 0.1 INPUT
1883	!
1884	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
1885	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
1886	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
1887	!
1888	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
1889	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
1890	!
1891	! 0.2 OUTPUT
1892	!
1893	REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) ! lai read variable and re-dimensioned
1894	!
1895	! 0.3 LOCAL
1896	!
1897	!
1898	CHARACTER(LEN=80) :: filename
1899	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
1900	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu
1901	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon
1902	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
1903	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
1904	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: laimap_lu
1905	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1906	!
1907	REAL(r_std) :: totarea
1908	INTEGER(i_std) :: idi, nbvmax
1909	CHARACTER(LEN=30) :: callsign
1910	!
1911	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
1912	!
1913	INTEGER :: ALLOC_ERR
1914	!
1915	!Config Key = LAI_FILE
1916	!Config Desc = Name of file from which the vegetation map is to be read
1917	!Config If = LAI_MAP
1918	!Config Def = ../surfmap/lai2D.nc
1919	!Config Help = The name of the file to be opened to read the LAI
1920	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
1921	!Config map which is derived from a Nicolas VIOVY one.
1922	!
1923	filename = 'lai2D.nc'
1924	CALL getin_p('LAI_FILE',filename)
1925	!
1926	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
1927	CALL bcast(iml)
1928	CALL bcast(jml)
1929	CALL bcast(lml)
1930	CALL bcast(tml)
1931	!
1932	!
1933	ALLOC_ERR=-1
1934	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
1935	IF (ALLOC_ERR/=0) THEN
1936	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
1937	STOP
1938	ENDIF
1939	ALLOC_ERR=-1
1940	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
1941	IF (ALLOC_ERR/=0) THEN
1942	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
1943	STOP
1944	ENDIF
1945	ALLOC_ERR=-1
1946	ALLOCATE(laimap_lu(iml,jml,nvm,tml), STAT=ALLOC_ERR)
1947	IF (ALLOC_ERR/=0) THEN
1948	WRITE(numout,*) "ERROR IN ALLOCATION of laimap_lu : ",ALLOC_ERR
1949	STOP
1950	ENDIF
1951	!
1952	!
1953	IF (is_root_prc) THEN
1954	CALL flinget(fid, 'longitude', iml, 0, 0, 0, 1, 1, lon_lu)
1955	CALL flinget(fid, 'latitude', jml, 0, 0, 0, 1, 1, lat_lu)
1956	CALL flinget(fid, 'LAI', iml, jml, nvm, tml, 1, 12, laimap_lu)
1957	!
1958	WHERE (laimap_lu(:,:,:,:) < zero )
1959	laimap_lu(:,:,:,:) = zero
1960	ENDWHERE
1961	!
1962	CALL flinclo(fid)
1963	ENDIF
1964	CALL bcast(lon_lu)
1965	CALL bcast(lat_lu)
1966	CALL bcast(laimap_lu)
1967	!
1968	ALLOC_ERR=-1
1969	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
1970	IF (ALLOC_ERR/=0) THEN
1971	WRITE(numout,*) "ERROR IN ALLOCATION of lon : ",ALLOC_ERR
1972	STOP
1973	ENDIF
1974	ALLOC_ERR=-1
1975	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
1976	IF (ALLOC_ERR/=0) THEN
1977	WRITE(numout,*) "ERROR IN ALLOCATION of lat : ",ALLOC_ERR
1978	STOP
1979	ENDIF
1980	!
1981	DO ip=1,iml
1982	lat(ip,:) = lat_lu(:)
1983	ENDDO
1984	DO jp=1,jml
1985	lon(:,jp) = lon_lu(:)
1986	ENDDO
1987	!
1988	ALLOC_ERR=-1
1989	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
1990	IF (ALLOC_ERR/=0) THEN
1991	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
1992	STOP
1993	ENDIF
1994	!
1995	! Consider all points a priori
1996	!
1997	!!$ mask(:,:) = 1
1998	mask(:,:) = 0
1999	!
2000	DO ip=1,iml
2001	DO jp=1,jml
2002	!
2003	! Exclude the points where there is never a LAI value. It is probably
2004	! an ocean point.
2005	!
2006	!!$ IF (ABS(SUM(laimap_lu(ip,jp,:,:))) < EPSILON(laimap_lu) ) THEN
2007	!!$ mask(ip,jp) = 0
2008	!!$ ENDIF
2009	!
2010	IF ( ANY(laimap_lu(ip,jp,:,:) < 20.) ) THEN
2011	mask(ip,jp) = 1
2012	ENDIF
2013	!
2014	ENDDO
2015	ENDDO
2016	!
2017	nbvmax = 20
2018	!
2019	callsign = 'LAI map'
2020	!
2021	ok_interpol = .FALSE.
2022	DO WHILE ( .NOT. ok_interpol )
2023	WRITE(numout,*) "Projection arrays for ",callsign," : "
2024	WRITE(numout,*) "nbvmax = ",nbvmax
2025	!
2026	ALLOC_ERR=-1
2027	ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
2028	IF (ALLOC_ERR/=0) THEN
2029	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
2030	STOP
2031	ENDIF
2032	sub_index(:,:,:)=0
2033	ALLOC_ERR=-1
2034	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
2035	IF (ALLOC_ERR/=0) THEN
2036	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
2037	STOP
2038	ENDIF
2039	sub_area(:,:)=zero
2040	!
2041	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
2042	& iml, jml, lon, lat, mask, callsign, &
2043	& nbvmax, sub_index, sub_area, ok_interpol)
2044
2045	!
2046	IF ( .NOT. ok_interpol ) THEN
2047	DEALLOCATE(sub_area)
2048	DEALLOCATE(sub_index)
2049	ENDIF
2050	!
2051	nbvmax = nbvmax * 2
2052	ENDDO
2053	!
2054	laimap(:,:,:) = zero
2055	!
2056	DO ib=1,nbpt
2057	idi = COUNT(sub_area(ib,:) > zero)
2058	IF ( idi > 0 ) THEN
2059	totarea = zero
2060	DO jj=1,idi
2061	ip = sub_index(ib,jj,1)
2062	jp = sub_index(ib,jj,2)
2063	DO jv=1,nvm
2064	DO it=1,12
2065	laimap(ib,jv,it) = laimap(ib,jv,it) + laimap_lu(ip,jp,jv,it)*sub_area(ib,jj)
2066	ENDDO
2067	ENDDO
2068	totarea = totarea + sub_area(ib,jj)
2069	ENDDO
2070	!
2071	! Normalize
2072	!
2073	laimap(ib,:,:) = laimap(ib,:,:)/totarea
2074
2075	!!$ IF ( ANY( laimap(ib,:,:) > 10 ) ) THEN
2076	!!$ WRITE(numout,*) "For point ",ib
2077	!!$ WRITE(numout,*) lalo(ib,1), lalo(ib,2)
2078	!!$ WRITE(numout,*) "For ib=",ib," we have LAI ",laimap(ib,:,1:idi)
2079	!!$ WRITE(numout,*) "Detail of projection :"
2080	!!$ WRITE(numout,*) sub_area(ib,1:idi)
2081	!!$ WRITE(numout,*) "Total for projection :" ,totarea
2082	!!$ ENDIF
2083	!
2084	ELSE
2085	WRITE(numout,*) 'On point ', ib, ' no points where found for interpolating the LAI map.'
2086	WRITE(numout,*) 'Location : ', lalo(ib,2), lalo(ib,1)
2087	DO jv=1,nvm
2088	laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux
2089	ENDDO
2090	WRITE(numout,*) 'Solved by putting the average LAI for the PFT all year long'
2091	ENDIF
2092	ENDDO
2093	!
2094	WRITE(numout,*) 'slowproc_interlai : Interpolation Done'
2095	!
2096	!
2097	!
2098	DEALLOCATE(lat_lu)
2099	DEALLOCATE(lon_lu)
2100	DEALLOCATE(lon)
2101	DEALLOCATE(lat)
2102	DEALLOCATE(laimap_lu)
2103	DEALLOCATE(mask)
2104	DEALLOCATE(sub_area)
2105	DEALLOCATE(sub_index)
2106	!
2107	RETURN
2108	!
2109	END SUBROUTINE slowproc_interlai_NEW
2110	!!
2111	!! Interpolate a vegetation map (by pft)
2112	!!
2113	!MM modif TAG 1.4 :
2114	! SUBROUTINE slowproc_update(nbpt, lalo, resolution, vegetmap, frac_nobiomap)
2115	!MM modif TAG 1.9.2 :
2116	! SUBROUTINE slowproc_update(nbpt, lalo, neighbours, resolution, contfrac, vegetmax, frac_nobio, veget_year, init)
2117	SUBROUTINE slowproc_update(nbpt, lalo, neighbours, resolution, contfrac, &
2118	& veget_last, frac_nobio_last, &
2119	& veget_next, frac_nobio_next, veget_year, init)
2120	!
2121	!
2122	!
2123	! 0.1 INPUT
2124	!
2125	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs
2126	! to be interpolated
2127	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo ! Vector of latitude and longitudes (beware of the order !)
2128	!MM modif TAG 1.4 : add grid variables for aggregate
2129	INTEGER(i_std), DIMENSION(nbpt,8), INTENT(in) :: neighbours ! Vector of neighbours for each grid point
2130	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
2131	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution ! The size in km of each grid-box in X and Y
2132	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of continent in the grid
2133	!
2134	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last ! old max vegetfrac
2135	REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(in) :: frac_nobio_last ! old fraction of the mesh which is
2136	! covered by ice, lakes, ...
2137	!
2138	INTEGER(i_std), INTENT(in) :: veget_year ! first year for landuse (0 == NO TIME AXIS)
2139	LOGICAL, OPTIONAL, INTENT(in) :: init ! initialisation : in case of dgvm, it forces update of all PFTs
2140	!
2141	! 0.2 OUTPUT
2142	!
2143	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(out) :: veget_next ! new max vegetfrac
2144	REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next ! new fraction of the mesh which is
2145	! covered by ice, lakes, ...
2146	!
2147	! 0.3 LOCAL
2148	!
2149	!
2150	CHARACTER(LEN=80) :: filename
2151	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, inobio, jv
2152	INTEGER(i_std) :: nb_coord,nb_var, nb_gat,nb_dim
2153	REAL(r_std) :: date, dt
2154	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: itau
2155	REAL(r_std), DIMENSION(1) :: time_counter
2156	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu
2157	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w, i_d_w
2158	LOGICAL :: exv, l_ex
2159	!
2160	!MM modif TAG 1.4 : suppression of all time axis reading and interpolation, replaced by each year 2D reading.
2161	! REAL(r_std), INTENT(inout) :: vegetmap(nbpt,nvm,293) ! vegetfrac read variable and re-dimensioned
2162	! REAL(r_std), INTENT(inout) :: frac_nobiomap(nbpt,nnobio,293) ! Fraction of the mesh which is covered by ice, lakes, ...
2163	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: vegmap ! last coord is time with only one value = 1
2164	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: vegmap_1 ! last coord is time with only one value = 1 (IF VEGET_YEAR == 0 , NO TIME AXIS)
2165	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
2166	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
2167	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
2168	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2169	!
2170	REAL(r_std) :: sumf, err, norm
2171	REAL(r_std) :: totarea
2172	INTEGER(i_std) :: idi, nbvmax
2173	CHARACTER(LEN=30) :: callsign
2174	!
2175	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
2176	!
2177	! for DGVM case :
2178	REAL(r_std) :: sum_veg ! sum of vegets
2179	REAL(r_std) :: sum_nobio ! sum of nobios
2180	REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets
2181	REAL(r_std) :: rapport ! (S-B) / (S-A)
2182	LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones)
2183	! e.g. in case of DGVM and not init (optional parameter)
2184	!
2185	LOGICAL, PARAMETER :: debug = .FALSE.
2186	!
2187	INTEGER :: ALLOC_ERR
2188	!
2189	!Config Key = VEGETATION_FILE
2190	!Config Desc = Name of file from which the vegetation map is to be read
2191	!Config If = LAND_USE
2192	!Config Def = pft_new.nc
2193	!Config Help = The name of the file to be opened to read a vegetation
2194	!Config map (in pft) is to be given here.
2195	!
2196	filename = 'pft_new.nc'
2197	CALL getin_p('VEGETATION_FILE',filename)
2198	!
2199	IF (is_root_prc) THEN
2200	IF (debug) THEN
2201	WRITE(numout,*) "Entering slowproc_update. Debug mode."
2202	WRITE (*,'(/," --> fliodmpf")')
2203	CALL fliodmpf (TRIM(filename))
2204	WRITE (*,'(/," --> flioopfd")')
2205	ENDIF
2206	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
2207	IF (debug) THEN
2208	WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord
2209	WRITE (*,'(" Number of variables in the file : ",I2)') nb_var
2210	WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat
2211	ENDIF
2212	ENDIF
2213	CALL bcast(nb_coord)
2214	CALL bcast(nb_var)
2215	CALL bcast(nb_gat)
2216	IF ( veget_year > 0 ) THEN
2217	IF (is_root_prc) &
2218	CALL flioinqv (fid,v_n="time_counter",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2219	CALL bcast(l_d_w)
2220	tml=l_d_w(1)
2221	WRITE(numout,*) "tml =",tml
2222	ENDIF
2223	IF (is_root_prc) &
2224	CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2225	CALL bcast(l_d_w)
2226	iml=l_d_w(1)
2227	WRITE(numout,*) "iml =",iml
2228
2229	IF (is_root_prc) &
2230	CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2231	CALL bcast(l_d_w)
2232	jml=l_d_w(1)
2233	WRITE(numout,*) "jml =",jml
2234
2235	IF (is_root_prc) &
2236	CALL flioinqv (fid,v_n="veget",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2237	CALL bcast(l_d_w)
2238	lml=l_d_w(1)
2239
2240	IF (lml /= nvm) &
2241	CALL ipslerr (3,'slowproc_update', &
2242	& 'Problem with vegetation file for Land Use','lml /= nvm', &
2243	& '(number of pft must be equal)')
2244	!
2245	ALLOC_ERR=-1
2246	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
2247	IF (ALLOC_ERR/=0) THEN
2248	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
2249	STOP
2250	ENDIF
2251	ALLOC_ERR=-1
2252	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
2253	IF (ALLOC_ERR/=0) THEN
2254	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
2255	STOP
2256	ENDIF
2257
2258	IF ( veget_year > 0 ) THEN
2259	IF (tml > 0) THEN
2260	ALLOC_ERR=-1
2261	ALLOCATE(itau(tml), STAT=ALLOC_ERR)
2262	IF (ALLOC_ERR/=0) THEN
2263	WRITE(numout,*) "ERROR IN ALLOCATION of itau : ",ALLOC_ERR
2264	STOP
2265	ENDIF
2266	ENDIF
2267	!
2268	IF (is_root_prc) THEN
2269	IF (tml > 0) THEN
2270	CALL fliogstc (fid, t_axis=itau,x_axis=lon_lu,y_axis=lat_lu)
2271	IF (veget_year <= tml) THEN
2272	CALL fliogetv (fid,"time_counter",time_counter, start=(/ veget_year /), count=(/ 1 /))
2273	WRITE(numout,*) "slowproc_update LAND_USE : the date read for vegetmax is ",time_counter
2274	ELSE
2275	CALL fliogetv (fid,"time_counter",time_counter, start=(/ tml /), count=(/ 1 /))
2276	WRITE(numout,*) "slowproc_update LAND_USE : You try to update vegetmax with a the date greater than in the file."
2277	WRITE(numout,*) " We will keep the last one :",time_counter
2278	ENDIF
2279	ELSE
2280	CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu)
2281	ENDIF
2282	ENDIF
2283	ENDIF
2284	IF (tml > 0) THEN
2285	CALL bcast(itau)
2286	ENDIF
2287	CALL bcast(lon_lu)
2288	CALL bcast(lat_lu)
2289	!
2290	ALLOC_ERR=-1
2291	ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR)
2292	IF (ALLOC_ERR/=0) THEN
2293	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
2294	STOP
2295	ENDIF
2296	ALLOC_ERR=-1
2297	ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR)
2298	IF (ALLOC_ERR/=0) THEN
2299	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
2300	STOP
2301	ENDIF
2302	!
2303	DO ip=1,iml
2304	lon_ful(ip,:)=lon_lu(ip)
2305	ENDDO
2306	DO jp=1,jml
2307	lat_ful(:,jp)=lat_lu(jp)
2308	ENDDO
2309	!
2310	!
2311	WRITE(numout,*) 'Reading the LAND USE vegetation file'
2312	!
2313	ALLOC_ERR=-1
2314	ALLOCATE(vegmap(iml,jml,nvm,1), STAT=ALLOC_ERR)
2315	IF (ALLOC_ERR/=0) THEN
2316	WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR
2317	STOP
2318	ENDIF
2319	IF ( veget_year == 0 ) THEN
2320	IF (is_root_prc) THEN
2321	ALLOC_ERR=-1
2322	ALLOCATE(vegmap_1(iml,jml,nvm), STAT=ALLOC_ERR)
2323	IF (ALLOC_ERR/=0) THEN
2324	WRITE(numout,*) "ERROR IN ALLOCATION of vegmap_1 : ",ALLOC_ERR
2325	STOP
2326	ENDIF
2327	ENDIF
2328	ENDIF
2329	!
2330	!!$ CALL flinopen &
2331	!!$ & (filename, .FALSE., iml, jml, lml, lon_ful, lat_ful, &
2332	!!$ & lev_ful, tml, itau, date, dt, fid)
2333	!=> FATAL ERROR FROM ROUTINE flinopen
2334	! --> No time axis found
2335
2336	!MM modif TAG 1.4 :
2337	! CALL flinget(fid, 'lon', iml, 0, 0, 0, 1, 1, lon_lu)
2338	! CALL flinget(fid, 'lat', jml, 0, 0, 0, 1, 1, lat_lu)
2339	! CALL flinget(fid, 'maxvegetfrac', iml, jml, nvm, tml, 1, 293, vegmap_lu)
2340	!FATAL ERROR FROM ROUTINE flinopen
2341	! --> No variable lon
2342	! We get only the right year
2343	!!$ CALL flinget(fid, 'maxvegetfrac', iml, jml, nvm, tml, veget_year, veget_year, vegmap)
2344	!!$ !
2345	!!$ CALL flinclo(fid)
2346
2347	IF (is_root_prc) &
2348	CALL flioinqv (fid,"maxvegetfrac",exv,nb_dims=nb_dim,len_dims=l_d_w,id_dims=i_d_w)
2349	CALL bcast(exv)
2350	CALL bcast(nb_dim)
2351	CALL bcast(l_d_w)
2352	CALL bcast(i_d_w)
2353
2354	IF (exv) THEN
2355	IF (debug) THEN
2356	WRITE (numout,'(" Number of dimensions : ",I2)') nb_dim
2357	WRITE (numout,'(" Dimensions :",/,5(1X,I7,:))') l_d_w(1:nb_dim)
2358	WRITE (numout,'(" Identifiers :",/,5(1X,I7,:))') i_d_w(1:nb_dim)
2359	ENDIF
2360	!
2361	IF ( veget_year > 0 ) THEN
2362	IF (is_root_prc) THEN
2363	IF (veget_year <= tml) THEN
2364	CALL fliogetv (fid,"maxvegetfrac", vegmap, start=(/ 1, 1, 1, veget_year /), count=(/ iml, jml, nvm, 1 /))
2365	ELSE
2366	CALL fliogetv (fid,"maxvegetfrac", vegmap, start=(/ 1, 1, 1, tml /), count=(/ iml, jml, nvm, 1 /))
2367	ENDIF
2368	ENDIF
2369	ELSE
2370	IF (is_root_prc) THEN
2371	CALL fliogetv (fid,"maxvegetfrac", vegmap_1, start=(/ 1, 1, 1 /), count=(/ iml, jml, nvm /))
2372	vegmap(:,:,:,1)=vegmap_1(:,:,:)
2373	DEALLOCATE(vegmap_1)
2374	ENDIF
2375	ENDIF
2376	CALL bcast(vegmap)
2377	IF (is_root_prc) CALL flioclo (fid)
2378	ELSE
2379	CALL ipslerr (3,'slowproc_update', &
2380	& 'Problem with vegetation file for Land Use.', &
2381	& "Variable maxvegetfrac doesn't exist.", &
2382	& '(verify your land use file.)')
2383	ENDIF
2384	!
2385	! Mask of permitted variables.
2386	!
2387	ALLOC_ERR=-1
2388	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2389	IF (ALLOC_ERR/=0) THEN
2390	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
2391	STOP
2392	ENDIF
2393	!
2394	mask(:,:) = 0
2395	DO ip=1,iml
2396	DO jp=1,jml
2397	sum_veg=SUM(vegmap(ip,jp,:,1))
2398	IF ( sum_veg .GE. min_sechiba .AND. sum_veg .LE. 1.-1.e-7) THEN
2399	mask(ip,jp) = 1
2400	IF (debug) THEN
2401	WRITE(numout,*) "update : SUM(vegmap(",ip,jp,")) = ",sum_veg
2402	ENDIF
2403	ELSEIF ( sum_veg .GT. 1.-1.e-7 .AND. sum_veg .LE. 2.) THEN
2404	! normalization
2405	vegmap(ip,jp,:,1) = vegmap(ip,jp,:,1) / sum_veg
2406	mask(ip,jp) = 1
2407	IF (debug) THEN
2408	WRITE(numout,*) "update : SUM(vegmap(",ip,jp,"))_c = ",SUM(vegmap(ip,jp,:,1))
2409	ENDIF
2410	ENDIF
2411	ENDDO
2412	ENDDO
2413	!
2414	!
2415	! The number of maximum vegetation map points in the GCM grid should
2416	! also be computed and not imposed here.
2417	!
2418	nbvmax = 200
2419	!
2420	callsign="Land Use Vegetation map"
2421	!
2422	ok_interpol = .FALSE.
2423	DO WHILE ( .NOT. ok_interpol )
2424	WRITE(numout,*) "Projection arrays for ",callsign," : "
2425	WRITE(numout,*) "nbvmax = ",nbvmax
2426	!
2427	ALLOC_ERR=-1
2428	ALLOCATE(sub_index(nbpt, nbvmax,2), STAT=ALLOC_ERR)
2429	IF (ALLOC_ERR/=0) THEN
2430	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
2431	STOP
2432	ENDIF
2433	sub_index(:,:,:)=0
2434
2435	ALLOC_ERR=-1
2436	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
2437	IF (ALLOC_ERR/=0) THEN
2438	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
2439	STOP
2440	ENDIF
2441	sub_area(:,:)=zero
2442	!
2443	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
2444	& iml, jml, lon_ful, lat_ful, mask, callsign, &
2445	& nbvmax, sub_index, sub_area, ok_interpol)
2446	!
2447	IF ( .NOT. ok_interpol ) THEN
2448	DEALLOCATE(sub_area)
2449	DEALLOCATE(sub_index)
2450	ENDIF
2451	!
2452	nbvmax = nbvmax * 2
2453	ENDDO
2454	!
2455	! Compute the logical for partial (only anthropic) PTFs update
2456	IF (PRESENT(init)) THEN
2457	partial_update = control%ok_dgvm .AND. .NOT. init
2458	ELSE
2459	partial_update = control%ok_dgvm
2460	ENDIF
2461	!
2462	IF ( .NOT. partial_update ) THEN
2463	!
2464	veget_next(:,:)=zero
2465	!
2466	DO ib = 1, nbpt
2467	idi=1
2468	sumf=zero
2469	DO WHILE ( sub_area(ib,idi) > zero )
2470	ip = sub_index(ib,idi,1)
2471	jp = sub_index(ib,idi,2)
2472	veget_next(ib,:) = veget_next(ib,:) + vegmap(ip,jp,:,1)*sub_area(ib,idi)
2473	sumf=sumf + sub_area(ib,idi)
2474	idi = idi +1
2475	ENDDO
2476	!!$ !
2477	!!$ ! Limit the smalest vegetation fraction to 0.5%
2478	!!$ !
2479	!!$ DO jv = 1, nvm
2480	!!$ IF ( veget_next(ib,jv) .LT. min_vegfrac ) THEN
2481	!!$ veget_next(ib,jv) = zero
2482	!!$ ENDIF
2483	!!$ ENDDO
2484	!
2485	! Normalize
2486	!
2487	IF (sumf > min_sechiba) THEN
2488	veget_next(ib,:) = veget_next(ib,:) / sumf
2489	ELSE
2490	WRITE(numout,*) "slowproc_update : No land point in the map for point ",&
2491	ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2492	CALL ipslerr (2,'slowproc_update', &
2493	& 'Problem with vegetation file for Land Use.', &
2494	& "No land point in the map for point", &
2495	& 'Keep old values. (verify your land use file.)')
2496	!!$ CALL slowproc_nearest (iml, lon_ful, lat_ful, &
2497	!!$ lalo(ib,2), lalo(ib,1), inear)
2498	IF (PRESENT(init)) THEN
2499	IF (init) THEN
2500	veget_next(ib,1) = un
2501	veget_next(ib,2:nvm) = zero
2502	ELSE
2503	veget_next(ib,:) = veget_last(ib,:)
2504	ENDIF
2505	ELSE
2506	veget_next(ib,:) = veget_last(ib,:)
2507	ENDIF
2508	ENDIF
2509	!
2510	IF (debug) THEN
2511	WRITE(numout,*) "SUM(veget_next(",ib,")) = ",SUM(veget_next(ib,:))
2512	ENDIF
2513	ENDDO
2514	ELSE
2515	DO ib = 1, nbpt
2516	! last veget for this point
2517	sum_veg=SUM(veget_last(ib,:))
2518	IF (debug) THEN
2519	WRITE(,) "SUM(veget_last(",ib,")) = ",sum_veg
2520	ENDIF
2521	!
2522	! If the DGVM is activated, only anthropiques PFT are utpdated,
2523	! other are copied
2524	veget_next(ib,:) = veget_last(ib,:)
2525	!
2526	! natural ones are initialized to zero.
2527	DO jv = 2, nvm
2528	! If the DGVM is activated, only anthropiques PFT are utpdated
2529	IF ( .NOT. natural(jv) ) THEN
2530	veget_next(ib,jv) = zero
2531	ENDIF
2532	ENDDO
2533	!
2534	idi=1
2535	sumf=zero
2536	DO WHILE ( sub_area(ib,idi) > zero )
2537	ip = sub_index(ib,idi,1)
2538	jp = sub_index(ib,idi,2)
2539	! If the DGVM is activated, only anthropic PFTs are utpdated
2540	DO jv = 2, nvm
2541	IF ( .NOT. natural(jv) ) THEN
2542	veget_next(ib,jv) = veget_next(ib,jv) + vegmap(ip,jp,jv,1)*sub_area(ib,idi)
2543	ENDIF
2544	ENDDO
2545	sumf=sumf + sub_area(ib,idi)
2546	idi = idi +1
2547	ENDDO
2548	!!$ !
2549	!!$ ! Limit the smalest vegetation fraction to 0.5%
2550	!!$ !
2551	!!$ DO jv = 2, nvm
2552	!!$ ! On anthropic and natural PFTs ?
2553	!!$ IF ( veget_next(ib,jv) .LT. min_vegfrac ) THEN
2554	!!$ veget_next(ib,jv) = zero
2555	!!$ ENDIF
2556	!!$ ENDDO
2557	!
2558	! Normalize
2559	!
2560	! Proposition de Pierre :
2561	! apres modification de la surface des PFTs anthropiques,
2562	! on doit conserver la proportion des PFTs naturels.
2563	! ie la somme des vegets est conservee
2564	! et PFT naturel / (somme des vegets - somme des vegets anthropiques)
2565	! est conservee.
2566	! Sum veget_next = old (sum veget_next Naturel) + (sum veget_next Anthropic)
2567	! = new (sum veget_next Naturel) + (sum veget_next Anthropic)
2568	! a / (S-A) = e / (S-B) ; b/(S-A) = f/(S-B)
2569	IF (sumf > min_sechiba) THEN
2570	sumvAnthro_old = zero
2571	sumvAnthro = zero
2572	DO jv = 2, nvm
2573	IF ( .NOT. natural(jv) ) THEN
2574	veget_next(ib,jv) = veget_next(ib,jv) / sumf
2575	sumvAnthro = sumvAnthro + veget_last(ib,jv)
2576	sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv)
2577	ENDIF
2578	ENDDO
2579	! conservation :
2580	rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old )
2581	DO jv = 1, nvm
2582	IF ( natural(jv) ) THEN
2583	veget_next(ib,jv) = veget_last(ib,jv) * rapport
2584	ENDIF
2585	ENDDO
2586	! test
2587	IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > 10*EPSILON(un) ) THEN
2588	WRITE(numout,*) "No conservation of sum of veget for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2589	WRITE(numout,*) "last sum of veget ",sum_veg," new sum of veget ",SUM(veget_next(ib,:))," error : ",&
2590	& SUM(veget_next(ib,:)) - sum_veg
2591	WRITE(numout,*) "Anthropic modifications : last ",sumvAnthro_old," new ",sumvAnthro
2592	CALL ipslerr (3,'slowproc_update', &
2593	& 'No conservation of sum of veget_next', &
2594	& "The sum of veget_next is different after reading Land Use map.", &
2595	& '(verify the dgvm case model.)')
2596	ENDIF
2597	ELSE
2598	WRITE(numout,*) "No land point in the map for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2599	! CALL ipslerr (3,'slowproc_update', &
2600	CALL ipslerr (2,'slowproc_update', &
2601	& 'Problem with vegetation file for Land Use.', &
2602	& "No land point in the map for point", &
2603	& '(verify your land use file.)')
2604	veget_next(ib,:) = veget_last(ib,:)
2605	ENDIF
2606
2607	ENDDO
2608	ENDIF
2609	!
2610	frac_nobio_next (:,:) = un
2611	!
2612	!MM
2613	! Work only for one nnobio !! (ie ice)
2614	DO inobio=1,nnobio
2615	DO jv=1,nvm
2616	!
2617	DO ib = 1, nbpt
2618	frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv)
2619	ENDDO
2620	ENDDO
2621	!
2622	ENDDO
2623	!
2624	DO ib = 1, nbpt
2625	sum_veg = SUM(veget_next(ib,:))
2626	sum_nobio = SUM(frac_nobio_next(ib,:))
2627	IF (sum_nobio < 0.) THEN
2628	frac_nobio_next(ib,:) = zero
2629	veget_next(ib,1) = veget_next(ib,1) - sum_nobio
2630	sum_veg = SUM(veget_next(ib,:))
2631	ENDIF
2632	sumf = sum_veg + sum_nobio
2633	IF (sumf > min_sechiba) THEN
2634	veget_next(ib,:) = veget_next(ib,:) / sumf
2635	frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf
2636	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2637	err=norm-un
2638	IF (debug) &
2639	WRITE(numout,*) "ib ",ib," SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un, sumf",err,sumf
2640	IF (abs(err) > -EPSILON(un)) THEN
2641	!MM 1.9.3
2642	! IF (abs(err) > 0.) THEN
2643	IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN
2644	frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err
2645	ELSE
2646	veget_next(ib,1) = veget_next(ib,1) - err
2647	ENDIF
2648	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2649	err=norm-un
2650	IF (debug) &
2651	WRITE(numout,*) "ib ",ib," SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un",err
2652	IF (abs(err) > EPSILON(un)) THEN
2653	!MM 1.9.3
2654	! IF (abs(err) > 0.) THEN
2655	WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2656	WRITE(numout,*) " err(sum-1.) = ",abs(err)
2657	CALL ipslerr (2,'slowproc_update', &
2658	& 'Problem with sum vegetation + sum fracnobio for Land Use.', &
2659	& "sum not equal to 1.", &
2660	& '(verify your land use file.)')
2661	ENDIF
2662	ENDIF
2663	ELSE
2664	WRITE(numout,*) "No vegetation nor frac_nobio for point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2665	WRITE(numout,*)"Replaced by bare_soil !! "
2666	veget_next(ib,1) = un
2667	veget_next(ib,2:nvm) = zero
2668	frac_nobio_next(ib,:) = zero
2669	!!$ CALL ipslerr (3,'slowproc_update', &
2670	!!$ & 'Problem with vegetation file for Land Use.', &
2671	!!$ & "No vegetation nor frac_nobio for point ", &
2672	!!$ & '(verify your land use file.)')
2673	ENDIF
2674	ENDDO
2675	!
2676	WRITE(numout,*) 'slowproc_update : Interpolation Done'
2677	!
2678	DEALLOCATE(vegmap)
2679	DEALLOCATE(lat_lu,lon_lu)
2680	DEALLOCATE(lat_ful,lon_ful)
2681	DEALLOCATE(mask)
2682	DEALLOCATE(sub_index,sub_area)
2683	!
2684	RETURN
2685	!
2686	END SUBROUTINE slowproc_update
2687
2688	!!
2689	!! Interpolate the IGBP vegetation map to the grid of the model
2690	!MM TAG 1.6 model !
2691	!!
2692	SUBROUTINE slowproc_interpol_OLD(nbpt, lalo, neighbours, resolution, veget, frac_nobio )
2693	!
2694	!
2695	!
2696	! 0.1 INPUT
2697	!
2698	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
2699	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
2700	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
2701	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
2702	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
2703	!
2704	! 0.2 OUTPUT
2705	!
2706	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
2707	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
2708	!
2709	! 0.3 LOCAL
2710	!
2711	INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes
2712	!
2713	!
2714	CHARACTER(LEN=80) :: filename
2715	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid
2716	REAL(r_std) :: lev(1), date, dt, coslat, pi
2717	INTEGER(i_std) :: itau(1)
2718	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
2719	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
2720	INTEGER, DIMENSION(nbpt,nolson) :: n_origveg
2721	INTEGER, DIMENSION(nbpt) :: n_found
2722	REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg
2723	REAL(r_std) :: vegcorr(nolson,nvm)
2724	REAL(r_std) :: nobiocorr(nolson,nnobio)
2725	CHARACTER(LEN=80) :: meter
2726	REAL(r_std) :: prog, sumf
2727	LOGICAL :: found
2728	INTEGER :: idi, ilast, ii, jv, inear, iprog
2729	REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
2730	!
2731	CALL get_vegcorr (nolson,vegcorr,nobiocorr)
2732	!
2733	!Config Key = VEGETATION_FILE
2734	!Config Desc = Name of file from which the vegetation map is to be read
2735	!Config If = !IMPOSE_VEG
2736	!Config Def = ../surfmap/carteveg5km.nc
2737	!Config Help = The name of the file to be opened to read the vegetation
2738	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
2739	!Config map which is derived from the IGBP one. We assume that we have
2740	!Config a classification in 87 types. This is Olson modified by Viovy.
2741	!
2742	filename = '../surfmap/carteveg5km.nc'
2743	CALL getin_p('VEGETATION_FILE',filename)
2744	!
2745	if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
2746	CALL bcast(iml)
2747	CALL bcast(jml)
2748	CALL bcast(lml)
2749	CALL bcast(tml)
2750	!
2751	!
2752	ALLOCATE(lat_ful(iml))
2753	ALLOCATE(lon_ful(iml))
2754	ALLOCATE(vegmap(iml))
2755	!
2756	WRITE(numout,*) 'Reading the vegetation file'
2757	!
2758	IF (is_root_prc) THEN
2759	CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
2760	CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
2761	CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
2762	!
2763	CALL flinclo(fid)
2764	ENDIF
2765
2766	CALL bcast(lon_ful)
2767	CALL bcast(lat_ful)
2768	CALL bcast(vegmap)
2769
2770	!
2771	IF (MAXVAL(vegmap) .LT. nolson) THEN
2772	WRITE(,) 'WARNING -- WARNING'
2773	WRITE(,) 'The vegetation map has to few vegetation types.'
2774	WRITE(,) 'If you are lucky it will work but please check'
2775	ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
2776	WRITE(,) 'More vegetation types in file than the code can'
2777	WRITE(,) 'deal with.: ', MAXVAL(vegmap), nolson
2778	STOP 'slowproc_interpol'
2779	ENDIF
2780	!
2781	ALLOCATE(lon_up(nbpt))
2782	ALLOCATE(lon_low(nbpt))
2783	ALLOCATE(lat_up(nbpt))
2784	ALLOCATE(lat_low(nbpt))
2785	!
2786	DO ib =1, nbpt
2787	!
2788	! We find the 4 limits of the grid-box. As we transform the resolution of the model
2789	! into longitudes and latitudes we do not have the problem of periodicity.
2790	! coslat is a help variable here !
2791	!
2792	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
2793	!
2794	lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
2795	lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
2796	!
2797	coslat = pi/180. * R_Earth
2798	!
2799	lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
2800	lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
2801	!
2802	!
2803	veget(ib,:) = zero
2804	frac_nobio (ib,:) = zero
2805	!
2806	ENDDO
2807	!
2808	! Get the limits of the integration domaine so that we can speed up the calculations
2809	!
2810	domaine_lon_min = MINVAL(lon_low)
2811	domaine_lon_max = MAXVAL(lon_up)
2812	domaine_lat_min = MINVAL(lat_low)
2813	domaine_lat_max = MAXVAL(lat_up)
2814	!
2815	!!$ WRITE(,) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max
2816	!!$ WRITE(,) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max
2817	!
2818	! Ensure that the fine grid covers the whole domain
2819	WHERE ( lon_ful(:) .LT. domaine_lon_min )
2820	lon_ful(:) = lon_ful(:) + 360.
2821	ENDWHERE
2822	!
2823	WHERE ( lon_ful(:) .GT. domaine_lon_max )
2824	lon_ful(:) = lon_ful(:) - 360.
2825	ENDWHERE
2826	!
2827	WRITE(numout,*) 'Interpolating the vegetation map :'
2828	WRITE(numout,'(2a40)')'0%--------------------------------------', &
2829	& '------------------------------------100%'
2830	!
2831	ilast = 1
2832	n_origveg(:,:) = 0
2833	!
2834	DO ip=1,iml
2835	!
2836	! Give a progress meter
2837	!
2838	! prog = ip/float(iml)*79.
2839	! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN
2840	! meter(NINT(prog)+1:NINT(prog)+1) = 'x'
2841	! WRITE(numout, advance="no", FMT='(a)') ACHAR(13)
2842	! WRITE(numout, advance="no", FMT='(a80)') meter
2843	! ENDIF
2844	iprog = NINT(float(ip)/float(iml)79.) - NINT(float(ip-1)/float(iml)79.)
2845	IF ( iprog .NE. 0 ) THEN
2846	WRITE(numout,'(a1,$)') 'x'
2847	ENDIF
2848	!
2849	! Only start looking for its place in the smaler grid if we are within the domaine
2850	! That should speed up things !
2851	!
2852	IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
2853	( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
2854	( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
2855	( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
2856	!
2857	! look for point on GCM grid which this point on fine grid belongs to.
2858	! First look at the point on the model grid where we arrived just before. There is
2859	! a good chace that neighbouring points on the fine grid fall into the same model
2860	! grid box.
2861	!
2862	!
2863	! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE
2864	! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED.
2865	!
2866	IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. &
2867	( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
2868	( lat_ful(ip) .GT. lat_low(ilast) ) .AND. &
2869	( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
2870	!
2871	! We were lucky
2872	!
2873	n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1
2874	!
2875	ELSE
2876	!
2877	! Otherwise, look everywhere.
2878	! Begin close to last grid point.
2879	!
2880	found = .FALSE.
2881	idi = 1
2882	!
2883	DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
2884	!
2885	! forward and backward
2886	!
2887	DO ii = 1,2
2888	!
2889	IF ( ii .EQ. 1 ) THEN
2890	ib = ilast - idi
2891	ELSE
2892	ib = ilast + idi
2893	ENDIF
2894	!
2895	IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
2896	IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. &
2897	( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
2898	( lat_ful(ip) .GT. lat_low(ib) ) .AND. &
2899	( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
2900	!
2901	n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1
2902	ilast = ib
2903	found = .TRUE.
2904	!
2905	ENDIF
2906	ENDIF
2907	!
2908	ENDDO
2909	!
2910	idi = idi + 1
2911	!
2912	ENDDO
2913	!
2914	ENDIF ! lucky/not lucky
2915	!
2916	ENDIF ! in the domain
2917	ENDDO
2918
2919	!
2920	! Now we know how many points of which Olson type from the fine grid fall
2921	! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
2922	!
2923
2924	!
2925	! determine number of points of the fine grid which fall into each box of the
2926	! coarse grid
2927	!
2928	DO ib = 1, nbpt
2929	n_found(ib) = SUM( n_origveg(ib,:) )
2930	ENDDO
2931
2932	!
2933	! determine fraction of Olson vegetation type in each box of the coarse grid
2934	!
2935	DO vid = 1, nolson
2936	WHERE ( n_found(:) .GT. 0 )
2937	frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std)
2938	ELSEWHERE
2939	frac_origveg(:,vid) = zero
2940	ENDWHERE
2941	ENDDO
2942
2943	!
2944	! now finally calculate coarse vegetation map
2945	! Find which model vegetation corresponds to each Olson type
2946	!
2947	DO vid = 1, nolson
2948	!
2949	DO jv = 1, nvm
2950	veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
2951	ENDDO
2952	!
2953	DO jv = 1, nnobio
2954	frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
2955	ENDDO
2956	!
2957	ENDDO
2958	!
2959	!
2960	WRITE(numout,*)
2961	WRITE(numout,*) 'Interpolation Done'
2962	!
2963	! Clean up the point of the map
2964	!
2965	DO ib = 1, nbpt
2966	!
2967	! Let us see if all points found something in the 5km map !
2968	!
2969	IF ( n_found(ib) .EQ. 0 ) THEN
2970	!
2971	! Now we need to handle some exceptions
2972	!
2973	IF ( lalo(ib,1) .LT. -56.0) THEN
2974	! Antartica
2975	frac_nobio(ib,:) = zero
2976	frac_nobio(ib,iice) = un
2977	veget(ib,:) = zero
2978	!
2979	ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
2980	! Artica
2981	frac_nobio(ib,:) = zero
2982	frac_nobio(ib,iice) = un
2983	veget(ib,:) = zero
2984	!
2985	ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
2986	! Greenland
2987	frac_nobio(ib,:) = zero
2988	frac_nobio(ib,iice) = un
2989	veget(ib,:) = zero
2990	!
2991	ELSE
2992	!
2993	WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box'
2994	WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib)
2995	WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib)
2996	!
2997	WRITE(numout,*) 'Looking for nearest point on the 5 km map'
2998	CALL slowproc_nearest (iml, lon_ful, lat_ful, &
2999	lalo(ib,2), lalo(ib,1), inear)
3000	WRITE(numout,*) 'Coordinates of the nearest point:', &
3001	lon_ful(inear),lat_ful(inear)
3002	!
3003	DO jv = 1, nvm
3004	veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3005	ENDDO
3006	!
3007	DO jv = 1, nnobio
3008	frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3009	ENDDO
3010	!
3011	ENDIF
3012	!
3013	ENDIF
3014	!
3015	!
3016	! Limit the smalest vegetation fraction to 0.5%
3017	!
3018	DO vid = 1, nvm
3019	IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3020	veget(ib,vid) = zero
3021	ENDIF
3022	ENDDO
3023	!
3024	sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3025	frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3026	veget(ib,:) = veget(ib,:)/sumf
3027	!
3028	!
3029	ENDDO
3030	!
3031	DEALLOCATE(lon_up)
3032	DEALLOCATE(lon_low)
3033	DEALLOCATE(lat_up)
3034	DEALLOCATE(lat_low)
3035	DEALLOCATE(lat_ful)
3036	DEALLOCATE(lon_ful)
3037	DEALLOCATE(vegmap)
3038	!
3039	RETURN
3040	!
3041	END SUBROUTINE slowproc_interpol_OLD
3042	!!
3043	!! Interpolate the IGBP vegetation map to the grid of the model
3044	!!
3045	SUBROUTINE slowproc_interpol_NEW(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
3046	!
3047	!
3048	!
3049	! 0.1 INPUT
3050	!
3051	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3052	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
3053	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
3054	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
3055	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
3056	REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
3057	!
3058	! 0.2 OUTPUT
3059	!
3060	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
3061	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
3062	!
3063	LOGICAL :: ok_interpol ! optionnal return of aggregate_vec
3064	!
3065	! 0.3 LOCAL
3066	!
3067	INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes
3068	!
3069	!
3070	CHARACTER(LEN=80) :: filename
3071	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid
3072	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
3073	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
3074	INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index
3075	REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg
3076	REAL(r_std), DIMENSION(nbpt) :: n_found
3077	REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg
3078	REAL(r_std) :: vegcorr(nolson,nvm)
3079	REAL(r_std) :: nobiocorr(nolson,nnobio)
3080	CHARACTER(LEN=40) :: callsign
3081	REAL(r_std) :: sumf, resol_lon, resol_lat
3082	INTEGER(i_std) :: idi, jv, inear, nbvmax
3083	!
3084	INTEGER :: ALLOC_ERR
3085	!
3086	n_origveg(:,:) = zero
3087	n_found(:) = zero
3088	!
3089	CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3090	!
3091	!Config Key = VEGETATION_FILE
3092	!Config Desc = Name of file from which the vegetation map is to be read
3093	!Config If = !IMPOSE_VEG
3094	!Config If = !LAND_USE
3095	!Config Def = ../surfmap/carteveg5km.nc
3096	!Config Help = The name of the file to be opened to read the vegetation
3097	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
3098	!Config map which is derived from the IGBP one. We assume that we have
3099	!Config a classification in 87 types. This is Olson modified by Viovy.
3100	!
3101	filename = '../surfmap/carteveg5km.nc'
3102	CALL getin_p('VEGETATION_FILE',filename)
3103	!
3104	if (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
3105	CALL bcast(iml)
3106	CALL bcast(jml)
3107	CALL bcast(lml)
3108	CALL bcast(tml)
3109	!
3110	!
3111	ALLOC_ERR=-1
3112	ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR)
3113	IF (ALLOC_ERR/=0) THEN
3114	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
3115	STOP
3116	ENDIF
3117	ALLOC_ERR=-1
3118	ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR)
3119	IF (ALLOC_ERR/=0) THEN
3120	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
3121	STOP
3122	ENDIF
3123	ALLOC_ERR=-1
3124	ALLOCATE(vegmap(iml), STAT=ALLOC_ERR)
3125	IF (ALLOC_ERR/=0) THEN
3126	WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR
3127	STOP
3128	ENDIF
3129	!
3130	WRITE(numout,*) 'Reading the OLSON type vegetation file'
3131	!
3132	IF (is_root_prc) THEN
3133	CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
3134	CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
3135	CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
3136	!
3137	CALL flinclo(fid)
3138	ENDIF
3139
3140	CALL bcast(lon_ful)
3141	CALL bcast(lat_ful)
3142	CALL bcast(vegmap)
3143
3144	!
3145	IF (MAXVAL(vegmap) .LT. nolson) THEN
3146	WRITE(numout,*) 'WARNING -- WARNING'
3147	WRITE(numout,*) 'The vegetation map has to few vegetation types.'
3148	WRITE(numout,*) 'If you are lucky it will work but please check'
3149	ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
3150	WRITE(numout,*) 'More vegetation types in file than the code can'
3151	WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson
3152	STOP 'slowproc_interpol'
3153	ENDIF
3154	!
3155	! Some assumptions on the vegetation file. This information should be
3156	! be computed or read from the file.
3157	! It is the reolution in meters of the grid of the vegetation file.
3158	!
3159	resol_lon = 5000.
3160	resol_lat = 5000.
3161	!
3162	! The number of maximum vegetation map points in the GCM grid should
3163	! also be computed and not imposed here.
3164	nbvmax = iml/nbpt
3165	!
3166	callsign="Vegetation map"
3167	!
3168	ok_interpol = .FALSE.
3169	DO WHILE ( .NOT. ok_interpol )
3170	WRITE(numout,*) "Projection arrays for ",callsign," : "
3171	WRITE(numout,*) "nbvmax = ",nbvmax
3172	!
3173	ALLOC_ERR=-1
3174	ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
3175	IF (ALLOC_ERR/=0) THEN
3176	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
3177	STOP
3178	ENDIF
3179	sub_index(:,:)=0
3180	ALLOC_ERR=-1
3181	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
3182	IF (ALLOC_ERR/=0) THEN
3183	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
3184	STOP
3185	ENDIF
3186	sub_area(:,:)=zero
3187	!
3188	CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, &
3189	& iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, &
3190	& nbvmax, sub_index, sub_area, ok_interpol)
3191	!
3192	IF ( .NOT. ok_interpol ) THEN
3193	DEALLOCATE(sub_area)
3194	DEALLOCATE(sub_index)
3195	!
3196	nbvmax = nbvmax * 2
3197	ELSE
3198	!
3199	DO ib = 1, nbpt
3200	idi=1
3201	DO WHILE ( sub_area(ib,idi) > zero )
3202	ip = sub_index(ib,idi)
3203	n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi)
3204	n_found(ib) = n_found(ib) + sub_area(ib,idi)
3205	idi = idi +1
3206	ENDDO
3207	ENDDO
3208	!
3209	ENDIF
3210	ENDDO
3211	!
3212	! Now we know how many points of which Olson type from the fine grid fall
3213	! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3214	!
3215	!
3216	! determine fraction of Olson vegetation type in each box of the coarse grid
3217	!
3218	DO vid = 1, nolson
3219	WHERE ( n_found(:) .GT. 0 )
3220	frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:)
3221	ELSEWHERE
3222	frac_origveg(:,vid) = zero
3223	ENDWHERE
3224	ENDDO
3225	!
3226	! now finally calculate coarse vegetation map
3227	! Find which model vegetation corresponds to each Olson type
3228	!
3229	veget(:,:) = zero
3230	frac_nobio(:,:) = zero
3231	!
3232	DO vid = 1, nolson
3233	!
3234	DO jv = 1, nvm
3235	veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
3236	ENDDO
3237	!
3238	DO jv = 1, nnobio
3239	frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
3240	ENDDO
3241	!
3242	ENDDO
3243	!
3244	WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
3245	!
3246	! Clean up the point of the map
3247	!
3248	DO ib = 1, nbpt
3249	!
3250	! Let us see if all points found something in the 5km map !
3251	!
3252	IF ( n_found(ib) .EQ. 0 ) THEN
3253	!
3254	! Now we need to handle some exceptions
3255	!
3256	IF ( lalo(ib,1) .LT. -56.0) THEN
3257	! Antartica
3258	frac_nobio(ib,:) = zero
3259	frac_nobio(ib,iice) = un
3260	veget(ib,:) = zero
3261	!
3262	ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3263	! Artica
3264	frac_nobio(ib,:) = zero
3265	frac_nobio(ib,iice) = un
3266	veget(ib,:) = zero
3267	!
3268	ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3269	! Greenland
3270	frac_nobio(ib,:) = zero
3271	frac_nobio(ib,iice) = un
3272	veget(ib,:) = zero
3273	!
3274	ELSE
3275	!
3276	WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib
3277	WRITE(numout,*) 'Longitude range : ', lalo(ib,2)
3278	WRITE(numout,*) 'Latitude range : ', lalo(ib,1)
3279	!
3280	WRITE(numout,*) 'Looking for nearest point on the 5 km map'
3281	CALL slowproc_nearest (iml, lon_ful, lat_ful, &
3282	lalo(ib,2), lalo(ib,1), inear)
3283	WRITE(numout,*) 'Coordinates of the nearest point:', &
3284	lon_ful(inear),lat_ful(inear)
3285	!
3286	DO jv = 1, nvm
3287	veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3288	ENDDO
3289	!
3290	DO jv = 1, nnobio
3291	frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3292	ENDDO
3293	!
3294	ENDIF
3295	!
3296	ENDIF
3297	!
3298	!
3299	! Limit the smalest vegetation fraction to 0.5%
3300	!
3301	DO vid = 1, nvm
3302	IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3303	veget(ib,vid) = zero
3304	ENDIF
3305	ENDDO
3306	!
3307	sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3308	frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3309	veget(ib,:) = veget(ib,:)/sumf
3310	!
3311	!
3312	ENDDO
3313	!
3314	DEALLOCATE(vegmap)
3315	DEALLOCATE(lat_ful, lon_ful)
3316	DEALLOCATE(sub_index)
3317	DEALLOCATE(sub_area)
3318
3319	!
3320	RETURN
3321	!
3322	END SUBROUTINE slowproc_interpol_NEW
3323
3324	!!
3325	!! Interpolate the IGBP vegetation map to the grid of the model
3326	!MM TAG 1.6 model !
3327	!!
3328	SUBROUTINE slowproc_interpol_OLD_g(nbpt, lalo, neighbours, resolution, veget, frac_nobio )
3329	!
3330	!
3331	!
3332	! 0.1 INPUT
3333	!
3334	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3335	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
3336	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
3337	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
3338	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
3339	!
3340	! 0.2 OUTPUT
3341	!
3342	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
3343	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
3344	!
3345	! 0.3 LOCAL
3346	!
3347	INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes
3348	!
3349	!
3350	CHARACTER(LEN=80) :: filename
3351	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, vid
3352	REAL(r_std) :: lev(1), date, dt, coslat
3353	INTEGER(i_std) :: itau(1)
3354	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
3355	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
3356	INTEGER, DIMENSION(nbpt,nolson) :: n_origveg
3357	INTEGER, DIMENSION(nbpt) :: n_found
3358	REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg
3359	REAL(r_std) :: vegcorr(nolson,nvm)
3360	REAL(r_std) :: nobiocorr(nolson,nnobio)
3361	CHARACTER(LEN=80) :: meter
3362	REAL(r_std) :: prog, sumf
3363	LOGICAL :: found
3364	INTEGER :: idi, ilast, ii, jv, inear, iprog
3365	REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
3366	!
3367	!
3368	CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3369	!
3370	!Config Key = VEGETATION_FILE
3371	!Config Desc = Name of file from which the vegetation map is to be read
3372	!Config If = !IMPOSE_VEG
3373	!Config Def = ../surfmap/carteveg5km.nc
3374	!Config Help = The name of the file to be opened to read the vegetation
3375	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
3376	!Config map which is derived from the IGBP one. We assume that we have
3377	!Config a classification in 87 types. This is Olson modified by Viovy.
3378	!
3379	filename = '../surfmap/carteveg5km.nc'
3380	CALL getin('VEGETATION_FILE',filename)
3381	!
3382	CALL flininfo(filename, iml, jml, lml, tml, fid)
3383	!
3384	!
3385	ALLOCATE(lat_ful(iml))
3386	ALLOCATE(lon_ful(iml))
3387	ALLOCATE(vegmap(iml))
3388	!
3389	WRITE(numout,*) 'Reading the vegetation file'
3390	!
3391	CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
3392	CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
3393	CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
3394	!
3395	CALL flinclo(fid)
3396
3397	!
3398	IF (MAXVAL(vegmap) .LT. nolson) THEN
3399	WRITE(,) 'WARNING -- WARNING'
3400	WRITE(,) 'The vegetation map has to few vegetation types.'
3401	WRITE(,) 'If you are lucky it will work but please check'
3402	ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
3403	WRITE(,) 'More vegetation types in file than the code can'
3404	WRITE(,) 'deal with.: ', MAXVAL(vegmap), nolson
3405	STOP 'slowproc_interpol'
3406	ENDIF
3407	!
3408	ALLOCATE(lon_up(nbpt))
3409	ALLOCATE(lon_low(nbpt))
3410	ALLOCATE(lat_up(nbpt))
3411	ALLOCATE(lat_low(nbpt))
3412	!
3413	DO ib =1, nbpt
3414	!
3415	! We find the 4 limits of the grid-box. As we transform the resolution of the model
3416	! into longitudes and latitudes we do not have the problem of periodicity.
3417	! coslat is a help variable here !
3418	!
3419	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
3420	!
3421	lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
3422	lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
3423	!
3424	coslat = pi/180. * R_Earth
3425	!
3426	lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
3427	lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
3428	!
3429	!
3430	veget(ib,:) = zero
3431	frac_nobio (ib,:) = zero
3432	!
3433	ENDDO
3434	!
3435	! Get the limits of the integration domaine so that we can speed up the calculations
3436	!
3437	domaine_lon_min = MINVAL(lon_low)
3438	domaine_lon_max = MAXVAL(lon_up)
3439	domaine_lat_min = MINVAL(lat_low)
3440	domaine_lat_max = MAXVAL(lat_up)
3441	!
3442	!!$ WRITE(,) 'DOMAINE lon :', domaine_lon_min, domaine_lon_max
3443	!!$ WRITE(,) 'DOMAINE lat :', domaine_lat_min, domaine_lat_max
3444	!
3445	! Ensure that the fine grid covers the whole domain
3446	WHERE ( lon_ful(:) .LT. domaine_lon_min )
3447	lon_ful(:) = lon_ful(:) + 360.
3448	ENDWHERE
3449	!
3450	WHERE ( lon_ful(:) .GT. domaine_lon_max )
3451	lon_ful(:) = lon_ful(:) - 360.
3452	ENDWHERE
3453	!
3454	WRITE(numout,*) 'Interpolating the vegetation map :'
3455	WRITE(numout,'(2a40)')'0%--------------------------------------', &
3456	& '------------------------------------100%'
3457	!
3458	ilast = 1
3459	n_origveg(:,:) = 0
3460	!
3461	DO ip=1,iml
3462	!
3463	! Give a progress meter
3464	!
3465	! prog = ip/float(iml)*79.
3466	! IF ( ABS(prog - NINT(prog)) .LT. 1/float(iml)*79. ) THEN
3467	! meter(NINT(prog)+1:NINT(prog)+1) = 'x'
3468	! WRITE(numout, advance="no", FMT='(a)') ACHAR(13)
3469	! WRITE(numout, advance="no", FMT='(a80)') meter
3470	! ENDIF
3471	iprog = NINT(float(ip)/float(iml)79.) - NINT(float(ip-1)/float(iml)79.)
3472	IF ( iprog .NE. 0 ) THEN
3473	WRITE(numout,'(a1,$)') 'x'
3474	ENDIF
3475	!
3476	! Only start looking for its place in the smaler grid if we are within the domaine
3477	! That should speed up things !
3478	!
3479	IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
3480	( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
3481	( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
3482	( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
3483	!
3484	! look for point on GCM grid which this point on fine grid belongs to.
3485	! First look at the point on the model grid where we arrived just before. There is
3486	! a good chace that neighbouring points on the fine grid fall into the same model
3487	! grid box.
3488	!
3489	!
3490	! THERE IS A BUG HERE !!! IF THE GCM GRID SITS ON THE DATE LINE WE WILL HAVE FOR INSTANCE
3491	! LON_LOW = -182 AND LON_UP = -178. THUS WE WILL ONLY PICK UP HALF THE POINTS NEEDED.
3492	!
3493	IF ( ( lon_ful(ip) .GT. lon_low(ilast) ) .AND. &
3494	( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
3495	( lat_ful(ip) .GT. lat_low(ilast) ) .AND. &
3496	( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
3497	!
3498	! We were lucky
3499	!
3500	n_origveg(ilast,NINT(vegmap(ip))) = n_origveg(ilast,NINT(vegmap(ip))) + 1
3501	!
3502	ELSE
3503	!
3504	! Otherwise, look everywhere.
3505	! Begin close to last grid point.
3506	!
3507	found = .FALSE.
3508	idi = 1
3509	!
3510	DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
3511	!
3512	! forward and backward
3513	!
3514	DO ii = 1,2
3515	!
3516	IF ( ii .EQ. 1 ) THEN
3517	ib = ilast - idi
3518	ELSE
3519	ib = ilast + idi
3520	ENDIF
3521	!
3522	IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
3523	IF ( ( lon_ful(ip) .GT. lon_low(ib) ) .AND. &
3524	( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
3525	( lat_ful(ip) .GT. lat_low(ib) ) .AND. &
3526	( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
3527	!
3528	n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + 1
3529	ilast = ib
3530	found = .TRUE.
3531	!
3532	ENDIF
3533	ENDIF
3534	!
3535	ENDDO
3536	!
3537	idi = idi + 1
3538	!
3539	ENDDO
3540	!
3541	ENDIF ! lucky/not lucky
3542	!
3543	ENDIF ! in the domain
3544	ENDDO
3545
3546	!
3547	! Now we know how many points of which Olson type from the fine grid fall
3548	! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3549	!
3550
3551	!
3552	! determine number of points of the fine grid which fall into each box of the
3553	! coarse grid
3554	!
3555	DO ib = 1, nbpt
3556	n_found(ib) = SUM( n_origveg(ib,:) )
3557	ENDDO
3558
3559	!
3560	! determine fraction of Olson vegetation type in each box of the coarse grid
3561	!
3562	DO vid = 1, nolson
3563	WHERE ( n_found(:) .GT. 0 )
3564	frac_origveg(:,vid) = REAL(n_origveg(:,vid),r_std) / REAL(n_found(:),r_std)
3565	ELSEWHERE
3566	frac_origveg(:,vid) = zero
3567	ENDWHERE
3568	ENDDO
3569
3570	!
3571	! now finally calculate coarse vegetation map
3572	! Find which model vegetation corresponds to each Olson type
3573	!
3574	DO vid = 1, nolson
3575	!
3576	DO jv = 1, nvm
3577	veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
3578	ENDDO
3579	!
3580	DO jv = 1, nnobio
3581	frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
3582	ENDDO
3583	!
3584	ENDDO
3585	!
3586	!
3587	WRITE(numout,*)
3588	WRITE(numout,*) 'Interpolation Done'
3589	!
3590	! Clean up the point of the map
3591	!
3592	DO ib = 1, nbpt
3593	!
3594	! Let us see if all points found something in the 5km map !
3595	!
3596	IF ( n_found(ib) .EQ. 0 ) THEN
3597	!
3598	! Now we need to handle some exceptions
3599	!
3600	IF ( lalo(ib,1) .LT. -56.0) THEN
3601	! Antartica
3602	frac_nobio(ib,:) = zero
3603	frac_nobio(ib,iice) = un
3604	veget(ib,:) = zero
3605	!
3606	ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3607	! Artica
3608	frac_nobio(ib,:) = zero
3609	frac_nobio(ib,iice) = un
3610	veget(ib,:) = zero
3611	!
3612	ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3613	! Greenland
3614	frac_nobio(ib,:) = zero
3615	frac_nobio(ib,iice) = un
3616	veget(ib,:) = zero
3617	!
3618	ELSE
3619	!
3620	WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box'
3621	WRITE(numout,*) 'Longitude range : ', lon_low(ib), lon_up(ib)
3622	WRITE(numout,*) 'Latitude range : ', lat_low(ib), lat_up(ib)
3623	!
3624	WRITE(numout,*) 'Looking for nearest point on the 5 km map'
3625	CALL slowproc_nearest (iml, lon_ful, lat_ful, &
3626	lalo(ib,2), lalo(ib,1), inear)
3627	WRITE(numout,*) 'Coordinates of the nearest point:', &
3628	lon_ful(inear),lat_ful(inear)
3629	!
3630	DO jv = 1, nvm
3631	veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3632	ENDDO
3633	!
3634	DO jv = 1, nnobio
3635	frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3636	ENDDO
3637	!
3638	ENDIF
3639	!
3640	ENDIF
3641	!
3642	!
3643	! Limit the smalest vegetation fraction to 0.5%
3644	!
3645	DO vid = 1, nvm
3646	IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3647	veget(ib,vid) = zero
3648	ENDIF
3649	ENDDO
3650	!
3651	sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3652	frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3653	veget(ib,:) = veget(ib,:)/sumf
3654	!
3655	!
3656	ENDDO
3657	!
3658	DEALLOCATE(lon_up)
3659	DEALLOCATE(lon_low)
3660	DEALLOCATE(lat_up)
3661	DEALLOCATE(lat_low)
3662	DEALLOCATE(lat_ful)
3663	DEALLOCATE(lon_ful)
3664	DEALLOCATE(vegmap)
3665	!
3666	RETURN
3667	!
3668	END SUBROUTINE slowproc_interpol_OLD_g
3669	!!
3670	!! Interpolate the IGBP vegetation map to the grid of the model
3671	!!
3672	SUBROUTINE slowproc_interpol_NEW_g(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
3673	!
3674	!
3675	!
3676	! 0.1 INPUT
3677	!
3678	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3679	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
3680	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
3681	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
3682	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
3683	REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
3684	!
3685	! 0.2 OUTPUT
3686	!
3687	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
3688	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
3689	!
3690	LOGICAL :: ok_interpol ! optionnal return of aggregate_vec
3691	!
3692	! 0.3 LOCAL
3693	!
3694	INTEGER(i_std), PARAMETER :: nolson = 94 ! Number of Olson classes
3695	!
3696	!
3697	CHARACTER(LEN=80) :: filename
3698	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid
3699	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
3700	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
3701	INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index
3702	REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg
3703	REAL(r_std), DIMENSION(nbpt) :: n_found
3704	REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg
3705	REAL(r_std) :: vegcorr(nolson,nvm)
3706	REAL(r_std) :: nobiocorr(nolson,nnobio)
3707	CHARACTER(LEN=40) :: callsign
3708	REAL(r_std) :: sumf, resol_lon, resol_lat
3709	INTEGER(i_std) :: idi, jv, inear, nbvmax
3710	!
3711	INTEGER :: ALLOC_ERR
3712	!
3713	n_origveg(:,:) = zero
3714	n_found(:) = zero
3715	!
3716	CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3717	!
3718	!Config Key = VEGETATION_FILE
3719	!Config Desc = Name of file from which the vegetation map is to be read
3720	!Config If = !IMPOSE_VEG
3721	!Config If = !LAND_USE
3722	!Config Def = ../surfmap/carteveg5km.nc
3723	!Config Help = The name of the file to be opened to read the vegetation
3724	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
3725	!Config map which is derived from the IGBP one. We assume that we have
3726	!Config a classification in 87 types. This is Olson modified by Viovy.
3727	!
3728	filename = '../surfmap/carteveg5km.nc'
3729	CALL getin('VEGETATION_FILE',filename)
3730	!
3731	CALL flininfo(filename, iml, jml, lml, tml, fid)
3732	!
3733	!
3734	ALLOC_ERR=-1
3735	ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR)
3736	IF (ALLOC_ERR/=0) THEN
3737	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
3738	STOP
3739	ENDIF
3740	ALLOC_ERR=-1
3741	ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR)
3742	IF (ALLOC_ERR/=0) THEN
3743	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
3744	STOP
3745	ENDIF
3746	ALLOC_ERR=-1
3747	ALLOCATE(vegmap(iml), STAT=ALLOC_ERR)
3748	IF (ALLOC_ERR/=0) THEN
3749	WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR
3750	STOP
3751	ENDIF
3752	!
3753	WRITE(numout,*) 'Reading the OLSON type vegetation file'
3754	!
3755	CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
3756	CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
3757	CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
3758	!
3759	CALL flinclo(fid)
3760	!
3761	IF (MAXVAL(vegmap) .LT. nolson) THEN
3762	WRITE(numout,*) 'WARNING -- WARNING'
3763	WRITE(numout,*) 'The vegetation map has to few vegetation types.'
3764	WRITE(numout,*) 'If you are lucky it will work but please check'
3765	ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
3766	WRITE(numout,*) 'More vegetation types in file than the code can'
3767	WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson
3768	STOP 'slowproc_interpol'
3769	ENDIF
3770	!
3771	! Some assumptions on the vegetation file. This information should be
3772	! be computed or read from the file.
3773	! It is the reolution in meters of the grid of the vegetation file.
3774	!
3775	resol_lon = 5000.
3776	resol_lat = 5000.
3777	!
3778	! The number of maximum vegetation map points in the GCM grid should
3779	! also be computed and not imposed here.
3780	nbvmax = iml/nbpt
3781	!
3782	callsign="Vegetation map"
3783	!
3784	ok_interpol = .FALSE.
3785	DO WHILE ( .NOT. ok_interpol )
3786	WRITE(numout,*) "Projection arrays for ",callsign," : "
3787	WRITE(numout,*) "nbvmax = ",nbvmax
3788	!
3789	ALLOC_ERR=-1
3790	ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
3791	IF (ALLOC_ERR/=0) THEN
3792	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
3793	STOP
3794	ENDIF
3795	sub_index(:,:)=0
3796	ALLOC_ERR=-1
3797	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
3798	IF (ALLOC_ERR/=0) THEN
3799	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
3800	STOP
3801	ENDIF
3802	sub_area(:,:)=zero
3803	!
3804	CALL aggregate (nbpt, lalo, neighbours, resolution, contfrac, &
3805	& iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, &
3806	& nbvmax, sub_index, sub_area, ok_interpol)
3807	!
3808	IF ( .NOT. ok_interpol ) THEN
3809	DEALLOCATE(sub_area)
3810	DEALLOCATE(sub_index)
3811	!
3812	nbvmax = nbvmax * 2
3813	ELSE
3814	!
3815	DO ib = 1, nbpt
3816	idi=1
3817	DO WHILE ( sub_area(ib,idi) > zero )
3818	ip = sub_index(ib,idi)
3819	n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi)
3820	n_found(ib) = n_found(ib) + sub_area(ib,idi)
3821	idi = idi +1
3822	ENDDO
3823	ENDDO
3824	!
3825	ENDIF
3826	ENDDO
3827	!
3828	! Now we know how many points of which Olson type from the fine grid fall
3829	! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3830	!
3831	!
3832	! determine fraction of Olson vegetation type in each box of the coarse grid
3833	!
3834	DO vid = 1, nolson
3835	WHERE ( n_found(:) .GT. 0 )
3836	frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:)
3837	ELSEWHERE
3838	frac_origveg(:,vid) = zero
3839	ENDWHERE
3840	ENDDO
3841	!
3842	! now finally calculate coarse vegetation map
3843	! Find which model vegetation corresponds to each Olson type
3844	!
3845	veget(:,:) = zero
3846	frac_nobio(:,:) = zero
3847	!
3848	DO vid = 1, nolson
3849	!
3850	DO jv = 1, nvm
3851	veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
3852	ENDDO
3853	!
3854	DO jv = 1, nnobio
3855	frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
3856	ENDDO
3857	!
3858	ENDDO
3859	!
3860	WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
3861	!
3862	! Clean up the point of the map
3863	!
3864	DO ib = 1, nbpt
3865	!
3866	! Let us see if all points found something in the 5km map !
3867	!
3868	IF ( n_found(ib) .EQ. 0 ) THEN
3869	!
3870	! Now we need to handle some exceptions
3871	!
3872	IF ( lalo(ib,1) .LT. -56.0) THEN
3873	! Antartica
3874	frac_nobio(ib,:) = zero
3875	frac_nobio(ib,iice) = un
3876	veget(ib,:) = zero
3877	!
3878	ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3879	! Artica
3880	frac_nobio(ib,:) = zero
3881	frac_nobio(ib,iice) = un
3882	veget(ib,:) = zero
3883	!
3884	ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3885	! Greenland
3886	frac_nobio(ib,:) = zero
3887	frac_nobio(ib,iice) = un
3888	veget(ib,:) = zero
3889	!
3890	ELSE
3891	!
3892	WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib
3893	WRITE(numout,*) 'Longitude range : ', lalo(ib,2)
3894	WRITE(numout,*) 'Latitude range : ', lalo(ib,1)
3895	!
3896	WRITE(numout,*) 'Looking for nearest point on the 5 km map'
3897	CALL slowproc_nearest (iml, lon_ful, lat_ful, &
3898	lalo(ib,2), lalo(ib,1), inear)
3899	WRITE(numout,*) 'Coordinates of the nearest point:', &
3900	lon_ful(inear),lat_ful(inear)
3901	!
3902	DO jv = 1, nvm
3903	veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3904	ENDDO
3905	!
3906	DO jv = 1, nnobio
3907	frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3908	ENDDO
3909	!
3910	ENDIF
3911	!
3912	ENDIF
3913	!
3914	!
3915	! Limit the smalest vegetation fraction to 0.5%
3916	!
3917	DO vid = 1, nvm
3918	IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3919	veget(ib,vid) = zero
3920	ENDIF
3921	ENDDO
3922	!
3923	sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3924	frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3925	veget(ib,:) = veget(ib,:)/sumf
3926	!
3927	!
3928	ENDDO
3929	!
3930	DEALLOCATE(vegmap)
3931	DEALLOCATE(lat_ful, lon_ful)
3932	DEALLOCATE(sub_index)
3933	DEALLOCATE(sub_area)
3934
3935	!
3936	RETURN
3937	!
3938	END SUBROUTINE slowproc_interpol_NEW_g
3939
3940
3941	!!
3942	!! looks for nearest grid point on the fine map
3943	!!
3944	SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
3945
3946	INTEGER(i_std), INTENT(in) :: iml
3947	REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5
3948	REAL(r_std), INTENT(in) :: lonmod, latmod
3949
3950	INTEGER(i_std), INTENT(out) :: inear
3951
3952	! REAL(r_std) :: pi
3953	REAL(r_std) :: pa, p
3954	REAL(r_std) :: coscolat, sincolat
3955	REAL(r_std) :: cospa, sinpa
3956	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
3957	INTEGER(i_std) :: i
3958	INTEGER(i_std), DIMENSION(1) :: ineartab
3959	INTEGER :: ALLOC_ERR
3960
3961	ALLOC_ERR=-1
3962	ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
3963	IF (ALLOC_ERR/=0) THEN
3964	WRITE(numout,*) "ERROR IN ALLOCATION of cosang : ",ALLOC_ERR
3965	STOP
3966	ENDIF
3967
3968	pa = pi/2.0 - latmod*pi/180.0 ! dist. entre pole n et point a
3969	cospa = COS(pa)
3970	sinpa = SIN(pa)
3971
3972	DO i = 1, iml
3973
3974	sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 )
3975	coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 )
3976
3977	p = (lonmod-lon5(i))*pi/180.0 ! angle entre a et b (leurs meridiens)
3978
3979	! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
3980	cosang(i) = cospacoscolat + sinpasincolat*COS(p)
3981
3982	ENDDO
3983
3984	ineartab = MAXLOC( cosang(:) )
3985	inear = ineartab(1)
3986
3987	DEALLOCATE(cosang)
3988	END SUBROUTINE slowproc_nearest
3989
3990	!!
3991	!! Interpolate the Zobler soil type map
3992	!!
3993	SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soiltype, clayfraction)
3994	!
3995	!
3996	! This subroutine should read the Zobler map and interpolate to the model grid. The method
3997	! is to get fraction of the three main soiltypes for each grid box.
3998	! The soil fraction are going to be put into the array soiltype in the following order :
3999	! coarse, medium and fine.
4000	!
4001	!
4002	! 0.1 INPUT
4003	!
4004	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
4005	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
4006	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
4007	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
4008	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
4009	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
4010	!
4011	! 0.2 OUTPUT
4012	!
4013	REAL(r_std), INTENT(out) :: soiltype(nbpt, nstm) ! Soil type map to be created from the Zobler map
4014	REAL(r_std), INTENT(out) :: clayfraction(nbpt) ! The fraction of clay as used by STOMATE
4015	!
4016	!
4017	! 0.3 LOCAL
4018	!
4019	INTEGER(i_std) :: nbvmax
4020	!
4021	CHARACTER(LEN=80) :: filename
4022	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, ilf, nbexp
4023	REAL(r_std) :: lev(1), date, dt
4024	INTEGER(i_std) :: itau(1)
4025	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, soiltext
4026	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
4027	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
4028	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
4029	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt
4030	REAL(r_std) :: sgn
4031	CHARACTER(LEN=30) :: callsign
4032	!
4033	! Number of texture classes in Zobler
4034	!
4035	INTEGER(i_std), PARAMETER :: classnb = 7
4036	REAL(r_std) :: textfrac_table(classnb, nstm)
4037	!
4038	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
4039	!
4040	INTEGER :: ALLOC_ERR
4041	!
4042	!
4043	CALL get_soilcorr (classnb, textfrac_table)
4044	!
4045	! Needs to be a configurable variable
4046	!
4047	!
4048	!Config Key = SOILTYPE_FILE
4049	!Config Desc = Name of file from which soil types are read
4050	!Config Def = ../surfmap/soils_param.nc
4051	!Config If = !IMPOSE_VEG
4052	!Config Help = The name of the file to be opened to read the soil types.
4053	!Config The data from this file is then interpolated to the grid of
4054	!Config of the model. The aim is to get fractions for sand loam and
4055	!Config clay in each grid box. This information is used for soil hydrology
4056	!Config and respiration.
4057	!
4058	filename = '../surfmap/soils_param.nc'
4059	CALL getin_p('SOILTYPE_FILE',filename)
4060	!
4061	IF (is_root_prc) THEN
4062	CALL flininfo(filename,iml, jml, lml, tml, fid)
4063	CALL flinclo(fid)
4064	ENDIF
4065	CALL bcast(iml)
4066	CALL bcast(jml)
4067	CALL bcast(lml)
4068	CALL bcast(tml)
4069	!
4070	! soils_param.nc file is 1° soit texture file.
4071	!
4072	ALLOC_ERR=-1
4073	ALLOCATE(lat_rel(iml,jml), STAT=ALLOC_ERR)
4074	IF (ALLOC_ERR/=0) THEN
4075	WRITE(numout,*) "ERROR IN ALLOCATION of lat_rel : ",ALLOC_ERR
4076	STOP
4077	ENDIF
4078	ALLOC_ERR=-1
4079	ALLOCATE(lon_rel(iml,jml), STAT=ALLOC_ERR)
4080	IF (ALLOC_ERR/=0) THEN
4081	WRITE(numout,*) "ERROR IN ALLOCATION of lon_rel : ",ALLOC_ERR
4082	STOP
4083	ENDIF
4084	ALLOC_ERR=-1
4085	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4086	IF (ALLOC_ERR/=0) THEN
4087	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
4088	STOP
4089	ENDIF
4090	ALLOC_ERR=-1
4091	ALLOCATE(soiltext(iml,jml), STAT=ALLOC_ERR)
4092	IF (ALLOC_ERR/=0) THEN
4093	WRITE(numout,*) "ERROR IN ALLOCATION of soiltext : ",ALLOC_ERR
4094	STOP
4095	ENDIF
4096	!
4097	IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
4098	CALL bcast(lon_rel)
4099	CALL bcast(lat_rel)
4100	CALL bcast(itau)
4101	CALL bcast(date)
4102	CALL bcast(dt)
4103
4104	!
4105	IF (is_root_prc) CALL flinget(fid, 'soiltext', iml, jml, lml, tml, 1, 1, soiltext)
4106	CALL bcast(soiltext)
4107	!
4108	IF (is_root_prc) CALL flinclo(fid)
4109	!
4110	nbexp = 0
4111	!
4112	!
4113	! Mask of permitted variables.
4114	!
4115	mask(:,:) = zero
4116	DO ip=1,iml
4117	DO jp=1,jml
4118	IF (soiltext(ip,jp) .GT. min_sechiba) THEN
4119	mask(ip,jp) = un
4120	ENDIF
4121	ENDDO
4122	ENDDO
4123	!
4124	nbvmax = 200
4125	!
4126	callsign = "Soil types"
4127	!
4128	ok_interpol = .FALSE.
4129	DO WHILE ( .NOT. ok_interpol )
4130	WRITE(numout,*) "Projection arrays for ",callsign," : "
4131	WRITE(numout,*) "nbvmax = ",nbvmax
4132	!
4133	ALLOC_ERR=-1
4134	ALLOCATE(solt(nbvmax), STAT=ALLOC_ERR)
4135	IF (ALLOC_ERR/=0) THEN
4136	WRITE(numout,*) "ERROR IN ALLOCATION of solt : ",ALLOC_ERR
4137	STOP
4138	ENDIF
4139	ALLOC_ERR=-1
4140	ALLOCATE(sub_index(nbpt,nbvmax,2), STAT=ALLOC_ERR)
4141	IF (ALLOC_ERR/=0) THEN
4142	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
4143	STOP
4144	ENDIF
4145	sub_index(:,:,:)=0
4146	ALLOC_ERR=-1
4147	ALLOCATE(sub_area(nbpt,nbvmax), STAT=ALLOC_ERR)
4148	IF (ALLOC_ERR/=0) THEN
4149	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
4150	STOP
4151	ENDIF
4152	sub_area(:,:)=zero
4153	!
4154	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
4155	& iml, jml, lon_rel, lat_rel, mask, callsign, &
4156	& nbvmax, sub_index, sub_area, ok_interpol)
4157	!
4158	IF ( .NOT. ok_interpol ) THEN
4159	DEALLOCATE(sub_area)
4160	DEALLOCATE(sub_index)
4161	DEALLOCATE(solt)
4162	!
4163	nbvmax = nbvmax * 2
4164	ENDIF
4165	ENDDO
4166	!
4167	DO ib =1, nbpt
4168	!
4169	soiltype(ib,:) = zero
4170	clayfraction(ib) = zero
4171	!
4172	! GO through the point we have found
4173	!
4174	!
4175	fopt = COUNT(sub_area(ib,:) > zero)
4176	!
4177	! Check that we found some points
4178	!
4179	IF ( fopt .EQ. 0) THEN
4180	nbexp = nbexp + 1
4181	soiltype(ib,:) = soiltype_default(:)
4182	clayfraction(ib) = clayfraction_default
4183	ELSE
4184	!
4185	DO ilf = 1,fopt
4186	solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
4187	ENDDO
4188	!
4189	sgn = zero
4190	!
4191	! Compute the average bare soil albedo parameters
4192	!
4193	DO ilf = 1,fopt
4194	!
4195	! We have to take care of two exceptions here : type 6 = glacier and type 0 = ocean
4196	!
4197	IF ( (solt(ilf) .LE. classnb) .AND. (solt(ilf) .GT. 0) .AND.&
4198	& (solt(ilf) .NE. 6)) THEN
4199	SELECTCASE(solt(ilf))
4200	CASE(1)
4201	soiltype(ib,1) = soiltype(ib,1) + sub_area(ib,ilf)
4202	CASE(2)
4203	soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4204	CASE(3)
4205	soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4206	CASE(4)
4207	soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4208	CASE(5)
4209	soiltype(ib,3) = soiltype(ib,3) + sub_area(ib,ilf)
4210	CASE(7)
4211	soiltype(ib,2) = soiltype(ib,2) + sub_area(ib,ilf)
4212	CASE DEFAULT
4213	WRITE(numout,*) 'We should not be here, an impossible case appeared'
4214	STOP 'slowproc_soilt'
4215	END SELECT
4216	clayfraction(ib) = clayfraction(ib) + &
4217	& textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
4218	sgn = sgn + sub_area(ib,ilf)
4219	ELSE
4220	IF (solt(ilf) .GT. classnb) THEN
4221	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
4222	STOP 'slowproc_soilt'
4223	ENDIF
4224	ENDIF
4225	!
4226	ENDDO
4227	!
4228	! Normalize the surface
4229	!
4230	IF ( sgn .LT. min_sechiba) THEN
4231	nbexp = nbexp + 1
4232	soiltype(ib,:) = soiltype_default(:)
4233	clayfraction(ib) = clayfraction_default
4234	ELSE
4235	soiltype(ib,:) = soiltype(ib,:)/sgn
4236	clayfraction(ib) = clayfraction(ib)/sgn
4237	ENDIF
4238	!
4239	ENDIF
4240	!
4241	ENDDO
4242	!
4243	IF ( nbexp .GT. 0 ) THEN
4244	WRITE(numout,*) 'slowproc_soilt : The interpolation of the bare soil albedo had ', nbexp
4245	WRITE(numout,*) 'slowproc_soilt : points without data. This are either coastal points or'
4246	WRITE(numout,*) 'slowproc_soilt : ice covered land.'
4247	WRITE(numout,*) 'slowproc_soilt : The problem was solved by using the default soil types.'
4248	ENDIF
4249	!
4250	DEALLOCATE (lat_rel)
4251	DEALLOCATE (lon_rel)
4252	DEALLOCATE (mask)
4253	DEALLOCATE (sub_area)
4254	DEALLOCATE (sub_index)
4255	DEALLOCATE (soiltext)
4256	DEALLOCATE (solt)
4257	!
4258	!
4259	RETURN
4260	!
4261	END SUBROUTINE slowproc_soilt
4262
4263
4264
4265	!! move from constantes_veg
4266
4267	!===
4268	SUBROUTINE get_vegcorr (nolson,vegcorr,nobiocorr)
4269	!---------------------------------------------------------------------
4270	INTEGER(i_std),INTENT(in) :: nolson
4271	REAL(r_std),DIMENSION(nolson,nvm),INTENT(out) :: vegcorr(nolson,nvm)
4272	REAL(r_std),DIMENSION(nolson,nnobio),INTENT(out) :: nobiocorr
4273	!-
4274	INTEGER(i_std) :: ib
4275	INTEGER(i_std),PARAMETER :: nolson94 = 94
4276	INTEGER(i_std),PARAMETER :: nvm13 = 13
4277	!---------------------------------------------------------------------
4278	IF (nolson /= nolson94) THEN
4279	WRITE(numout,*) nolson,nolson94
4280	CALL ipslerr(3,'get_vegcorr', '', '',&
4281	& 'wrong number of OLSON vegetation types.')
4282	ENDIF
4283	IF (nvm /= nvm13) THEN
4284	WRITE(numout,*) nvm,nvm13
4285	CALL ipslerr(3,'get_vegcorr', '', '',&
4286	& 'wrong number of SECHIBA vegetation types.')
4287	ENDIF
4288	!-
4289	! 1 set the indices of non-biospheric surface types to 0.
4290	!-
4291	nobiocorr(:,:) = 0.
4292	!-
4293	! 2 Here we construct the correspondance table
4294	! between Olson and the following SECHIBA Classes.
4295	! vegcorr(i,:)+nobiocorr(i,:) = 1. for all i.
4296	!-
4297	! The modified OLSON types found in file carteveg5km.nc
4298	! created by Nicolas Viovy :
4299	! 1 Urban
4300	vegcorr( 1,:) = &
4301	& (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4302	! 2 Cool low sparse grassland
4303	vegcorr( 2,:) = &
4304	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4305	! 3 Cold conifer forest
4306	vegcorr( 3,:) = &
4307	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4308	! 4 Cold deciduous conifer forest
4309	vegcorr( 4,:) = &
4310	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0/)
4311	! 5 Cool Deciduous broadleaf forest
4312	vegcorr( 5,:) = &
4313	& (/0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4314	! 6 Cool evergreen broadleaf forests
4315	vegcorr( 6,:) = &
4316	& (/0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4317	! 7 Cool tall grasses and shrubs
4318	vegcorr( 7,:) = &
4319	& (/0.1, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4320	! 8 Warm C3 tall grasses and shrubs
4321	vegcorr( 8,:) = &
4322	& (/0.1, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4323	! 9 Warm C4 tall grases and shrubs
4324	vegcorr( 9,:) = &
4325	& (/0.1, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0/)
4326	! 10 Bare desert
4327	vegcorr(10,:) = &
4328	& (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4329	! 11 Cold upland tundra
4330	vegcorr(11,:) = &
4331	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4332	! 12 Cool irrigated grassland
4333	vegcorr(12,:) = &
4334	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0, 0.0/)
4335	! 13 Semi desert
4336	vegcorr(13,:) = &
4337	& (/0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0/)
4338	! 14 Glacier ice
4339	vegcorr(14,:) = &
4340	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4341	nobiocorr(14,iice) = 1.
4342	! 15 Warm wooded wet swamp
4343	vegcorr(15,:) = &
4344	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0/)
4345	! 16 Inland water
4346	vegcorr(16,:) = &
4347	& (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4348	! 17 sea water
4349	vegcorr(17,:) = &
4350	& (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4351	! 18 cool shrub evergreen
4352	vegcorr(18,:) = &
4353	& (/0.1, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0, 0.0/)
4354	! 19 cold shrub deciduous
4355	vegcorr(19,:) = &
4356	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.6, 0.0, 0.0, 0.0/)
4357	! 20 Cold evergreen forest and fields
4358	vegcorr(20,:) = &
4359	& (/0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0/)
4360	! 21 cool rain forest
4361	vegcorr(21,:) = &
4362	& (/0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4363	! 22 cold conifer boreal forest
4364	vegcorr(22,:) = &
4365	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4366	! 23 cool conifer forest
4367	vegcorr(23,:) = &
4368	& (/0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4369	! 24 warm mixed forest
4370	vegcorr(24,:) = &
4371	& (/0.0, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0/)
4372	! 25 cool mixed forest
4373	vegcorr(25,:) = &
4374	& (/0.0, 0.0, 0.0, 0.4, 0.0, 0.4, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4375	! 26 cool broadleaf forest
4376	vegcorr(26,:) = &
4377	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0/)
4378	! 27 cool deciduous broadleaf forest
4379	vegcorr(27,:) = &
4380	& (/0.0, 0.0, 0.0, 0.0, 0.3, 0.5, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4381	! 28 warm montane tropical forest
4382	vegcorr(28,:) = &
4383	& (/0.0, 0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0/)
4384	! 29 warm seasonal tropical forest
4385	vegcorr(29,:) = &
4386	& (/0.0, 0.5, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0/)
4387	! 30 cool crops and towns
4388	vegcorr(30,:) = &
4389	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0/)
4390	! 31 warm crops and towns
4391	vegcorr(31,:) = &
4392	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8/)
4393	! 32 cool crops and towns
4394	vegcorr(32,:) = &
4395	& (/0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0/)
4396	! 33 warm dry tropical woods
4397	vegcorr(33,:) = &
4398	& (/0.2, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0/)
4399	! 34 warm tropical rain forest
4400	vegcorr(34,:) = &
4401	& (/0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4402	! 35 warm tropical degraded forest
4403	vegcorr(35,:) = &
4404	& (/0.1, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0/)
4405	! 36 warm corn and beans cropland
4406	vegcorr(36,:) = &
4407	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/)
4408	! 37 cool corn and bean cropland
4409	vegcorr(37,:) = &
4410	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0/)
4411	! 38 warm rice paddy and field
4412	vegcorr(38,:) = &
4413	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/)
4414	! 39 hot irrigated cropland
4415	vegcorr(39,:) = &
4416	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/)
4417	! 40 cool irrigated cropland
4418	vegcorr(40,:) = &
4419	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0/)
4420	! 41 cold irrigated cropland
4421	vegcorr(41,:) = &
4422	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0/)
4423	! 42 cool grasses and shrubs
4424	vegcorr(42,:) = &
4425	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0, 0.7, 0.0, 0.0, 0.0/)
4426	! 43 hot and mild grasses and shrubs
4427	vegcorr(43,:) = &
4428	& (/0.2, 0.0, 0.1, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0/)
4429	! 44 cold grassland
4430	vegcorr(44,:) = &
4431	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0, 0.0/)
4432	! 45 Savanna (woods) C3
4433	vegcorr(45,:) = &
4434	& (/0.1, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7, 0.0, 0.0, 0.0/)
4435	! 46 Savanna woods C4
4436	vegcorr(46,:) = &
4437	& (/0.1, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7, 0.0, 0.0/)
4438	! 47 Mire, bog, fen
4439	vegcorr(47,:) = &
4440	& (/0.1, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7, 0.0, 0.0, 0.0/)
4441	! 48 Warm marsh wetland
4442	vegcorr(48,:) = &
4443	& (/0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4444	! 49 cold marsh wetland
4445	vegcorr(49,:) = &
4446	& (/0.0, 0.0, 0.0, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4447	! 50 mediteraean scrub
4448	vegcorr(50,:) = &
4449	& (/0.1, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4450	! 51 Cool dry woody scrub
4451	vegcorr(51,:) = &
4452	& (/0.3, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0, 0.0/)
4453	! 52 Warm dry evergreen woods
4454	vegcorr(52,:) = &
4455	& (/0.1, 0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4456	! 53 Volcanic rocks
4457	vegcorr(53,:) = &
4458	& (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4459	! 54 sand desert
4460	vegcorr(54,:) = &
4461	& (/1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4462	! 55 warm semi desert shrubs
4463	vegcorr(55,:) = &
4464	& (/0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0/)
4465	! 56 cool semi desert shrubs
4466	vegcorr(56,:) = &
4467	& (/0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0/)
4468	! 57 semi desert sage
4469	vegcorr(57,:) = &
4470	& (/0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0, 0.0/)
4471	! 58 Barren tundra
4472	vegcorr(58,:) = &
4473	& (/0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0/)
4474	! 59 cool southern hemisphere mixed forest
4475	vegcorr(59,:) = &
4476	& (/0.1, 0.0, 0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0/)
4477	! 60 cool fields and woods
4478	vegcorr(60,:) = &
4479	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0/)
4480	! 61 warm forest and filed
4481	vegcorr(61,:) = &
4482	& (/0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6/)
4483	! 62 cool forest and field
4484	vegcorr(62,:) = &
4485	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0/)
4486	! 63 warm C3 fields and woody savanna
4487	vegcorr(63,:) = &
4488	& (/0.1, 0.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0/)
4489	! 64 warm C4 fields and woody savanna
4490	vegcorr(64,:) = &
4491	& (/0.1, 0.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6/)
4492	! 65 cool fields and woody savanna
4493	vegcorr(65,:) = &
4494	& (/0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0/)
4495	! 66 warm succulent and thorn scrub
4496	vegcorr(66,:) = &
4497	& (/0.1, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0/)
4498	! 67 cold small leaf mixed woods
4499	vegcorr(67,:) = &
4500	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.3, 0.0, 0.5, 0.0, 0.0, 0.0/)
4501	! 68 cold deciduous and mixed boreal fores
4502	vegcorr(68,:) = &
4503	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0/)
4504	! 69 cold narrow conifers
4505	vegcorr(69,:) = &
4506	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0/)
4507	! 70 cold wooded tundra
4508	vegcorr(70,:) = &
4509	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.7, 0.0, 0.0, 0.0/)
4510	! 71 cold heath scrub
4511	vegcorr(71,:) = &
4512	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.7, 0.0, 0.0, 0.0/)
4513	! 72 Polar and alpine desert
4514	vegcorr(72,:) = &
4515	& (/0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0/)
4516	! 73 warm Mangrove
4517	vegcorr(73,:) = &
4518	& (/0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4519	! 74 cool crop and water mixtures
4520	vegcorr(74,:) = &
4521	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0/)
4522	! 75 cool southern hemisphere mixed forest
4523	vegcorr(75,:) = &
4524	& (/0.0, 0.0, 0.0, 0.0, 0.4, 0.4, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4525	! 76 cool moist eucalyptus
4526	vegcorr(76,:) = &
4527	& (/0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0/)
4528	! 77 warm rain green tropical forest
4529	vegcorr(77,:) = &
4530	& (/0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4531	! 78 warm C3 woody savanna
4532	vegcorr(78,:) = &
4533	& (/0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0, 0.0/)
4534	! 79 warm C4 woody savanna
4535	vegcorr(79,:) = &
4536	& (/0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0/)
4537	! 80 cool woody savanna
4538	vegcorr(80,:) = &
4539	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.0, 0.6, 0.0, 0.0, 0.0/)
4540	! 81 cold woody savanna
4541	vegcorr(81,:) = &
4542	& (/0.0, 0.0, 0.0, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6, 0.0, 0.0, 0.0/)
4543	! 82 warm broadleaf crops
4544	vegcorr(82,:) = &
4545	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0/)
4546	! 83 warm C3 grass crops
4547	vegcorr(83,:) = &
4548	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9, 0.0/)
4549	! 84 warm C4 grass crops
4550	vegcorr(84,:) = &
4551	& (/0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9/)
4552	! 85 cool grass crops
4553	vegcorr(85,:) = &
4554	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0/)
4555	! 86 warm C3 crops grass,shrubs
4556	vegcorr(86,:) = &
4557	& (/0.0, 0.0, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0/)
4558	! 87 cool crops,grass,shrubs
4559	vegcorr(87,:) = &
4560	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0/)
4561	! 88 warm evergreen tree crop
4562	vegcorr(88,:) = &
4563	& (/0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2/)
4564	! 89 cool evergreen tree crop
4565	vegcorr(89,:) = &
4566	& (/0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0/)
4567	! 90 cold evergreen tree crop
4568	vegcorr(90,:) = &
4569	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0/)
4570	! 91 warm deciduous tree crop
4571	vegcorr(91,:) = &
4572	& (/0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2/)
4573	! 92 cool deciduous tree crop
4574	vegcorr(92,:) = &
4575	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.2, 0.0/)
4576	! 93 cold deciduous tree crop
4577	vegcorr(93,:) = &
4578	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8, 0.0, 0.0, 0.0, 0.2, 0.0/)
4579	! 94 wet sclerophylic forest
4580	vegcorr(94,:) = &
4581	& (/0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4582	!-
4583	! 3 Check the mapping for the Olson types which are going into the
4584	! the veget and nobio array.
4585	!-
4586	DO ib=1,nolson
4587	IF ( ABS(SUM(vegcorr(ib,:))+SUM(nobiocorr(ib,:))-1.0) &
4588	& > EPSILON(1.0)) THEN
4589	WRITE(numout,*) 'Wrong correspondance for Olson type :', ib
4590	CALL ipslerr(3,'get_vegcorr', '', '',&
4591	& 'Wrong correspondance for Olson type.')
4592	ENDIF
4593	ENDDO
4594	!-------------------------
4595	END SUBROUTINE get_vegcorr
4596
4597
4598
4599	!
4600	SUBROUTINE get_soilcorr (nzobler,textfrac_table)
4601	!!--------------------------------------------------------------------
4602	!! The "get_soilcorr" routine defines the table of correspondence
4603	!! between the Zobler types and the three texture
4604	!! types known by SECHIBA & STOMATE : silt, sand and clay
4605	!!--------------------------------------------------------------------
4606	INTEGER(i_std),INTENT(in) :: nzobler
4607	REAL(r_std),DIMENSION(nzobler,nstm),INTENT(out) :: textfrac_table
4608	!-
4609	INTEGER(i_std),PARAMETER :: nbtypes_zobler = 7
4610	INTEGER(i_std) :: ib
4611	!---------------------------------------------------------------------
4612	IF (nzobler /= nbtypes_zobler) THEN
4613	CALL ipslerr(3,'get_soilcorr', 'nzobler /= nbtypes_zobler',&
4614	& 'We do not have the correct number of classes', &
4615	& ' in the code for the file.')
4616	ENDIF
4617	!-
4618	! Textural fraction for : silt sand clay
4619	!-
4620	textfrac_table(1,:) = (/ 0.12, 0.82, 0.06 /)
4621	textfrac_table(2,:) = (/ 0.32, 0.58, 0.10 /)
4622	textfrac_table(3,:) = (/ 0.39, 0.43, 0.18 /)
4623	textfrac_table(4,:) = (/ 0.15, 0.58, 0.27 /)
4624	textfrac_table(5,:) = (/ 0.34, 0.32, 0.34 /)
4625	textfrac_table(6,:) = (/ 0.00, 1.00, 0.00 /)
4626	textfrac_table(7,:) = (/ 0.39, 0.43, 0.18 /)
4627
4628	DO ib=1,nzobler
4629	IF (ABS(SUM(textfrac_table(ib,:))-1.0) > EPSILON(1.0)) THEN
4630	WRITE(numout,*) &
4631	& 'Error in the correspondence table', &
4632	& ' sum is not equal to 1 in', ib
4633	WRITE(numout,*) textfrac_table(ib,:)
4634	CALL ipslerr(3,'get_soilcorr', 'SUM(textfrac_table(ib,:)) /= 1.0',&
4635	& '', 'Error in the correspondence table')
4636	ENDIF
4637	ENDDO
4638	!--------------------------
4639	END SUBROUTINE get_soilcorr
4640
4641
4642	!===
4643	FUNCTION tempfunc (temp_in) RESULT (tempfunc_result)
4644	!!--------------------------------------------------------------------
4645	!! FUNCTION tempfunc (temp_in) RESULT (tempfunc_result)
4646	!! this function interpolates value between ztempmin and ztempmax
4647	!! used for lai detection
4648	!!--------------------------------------------------------------------
4649	REAL(r_std),INTENT(in) :: temp_in !! Temperature in degre Kelvin
4650	REAL(r_std) :: tempfunc_result
4651	!-
4652	REAL(r_std),PARAMETER :: ztempmin=273._r_std !! Temperature for laimin
4653	REAL(r_std),PARAMETER :: ztempmax=293._r_std !! Temperature for laimax
4654	REAL(r_std) :: zfacteur !! Interpolation factor
4655	!---------------------------------------------------------------------
4656	zfacteur = un/(ztempmax-ztempmin)**2
4657	IF (temp_in > ztempmax) THEN
4658	tempfunc_result = un
4659	ELSEIF (temp_in < ztempmin) THEN
4660	tempfunc_result = zero
4661	ELSE
4662	tempfunc_result = un-zfacteur(ztempmax-temp_in)*2
4663	ENDIF
4664	!--------------------
4665	END FUNCTION tempfunc
4666
4667
4668	END MODULE slowproc

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

Download in other formats: