Context Navigation

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

Last change on this file since 311 was 311, checked in by didier.solyga, 13 years ago
Move and clean the rest of the externalized parameters from sechiba and stomate to src_parameters. Add two subroutines in constantes. Correct Olson type number 79 in vegcorr
File size: 169.8 KB

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

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

Download in other formats: