Context Navigation

source: perso/abdelouhab.djerrah/ORCHIDEE_ORIG/src_sechiba/slowproc.f90 @ 938

Last change on this file since 938 was 372, checked in by martial.mancip, 13 years ago
Algorithm in slowproc_update LAND USE + DGVM modified in case of reforestation.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 158.2 KB

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

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

Download in other formats: