Context Navigation

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

Last change on this file since 103 was 64, checked in by didier.solyga, 14 years ago
Import first version of ORCHIDEE_EXT
File size: 174.6 KB

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

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

Download in other formats: