Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/slowproc.f90 @ 6268

Last change on this file since 6268 was 846, checked in by didier.solyga, 12 years ago
Formatted labels so a script can automatically generate the orchidee.default file.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 176.1 KB

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

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

Download in other formats: