Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

slowproc.f90 @ 3850

Last change on this file since 3850 was 8, checked in by orchidee, 14 years ago
import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 156.6 KB

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

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

Context Navigation

source: tags/ORCHIDEE_1_9_5/ORCHIDEE/src_sechiba/slowproc.f90 @ 3850

Download in other formats: