Context Navigation

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

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

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

Download in other formats: