Context Navigation

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

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

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

Download in other formats: