Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE/src_sechiba/slowproc.f90

Last change on this file was 45, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

There is no change in numerical result after these commits.

Property svn:keywords set to HeadURL Date Author Revision

File size: 157.6 KB

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

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

Download in other formats: