Context Navigation

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

Last change on this file since 407 was 386, checked in by didier.solyga, 13 years ago
Correct2 wrong loops in slowproc when dgvm is activated. Replace PRINT instructions by WRITE instructions. Add a call to ipslerr in stomate alloc in case of wrong values of L0 and R0
File size: 169.9 KB

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

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

Download in other formats: