Context Navigation

slowproc.f90 @ 854

Last change on this file since 854 was 821, checked in by sebastiaan.luyssaert, 12 years ago
Bug fixed, reintroduced loop of PFTs that god lost during documentation effort
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 186.0 KB

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

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

Context Navigation

source: perso/abdelouhab.djerrah/ORCHIDEE/src_sechiba/slowproc.f90 @ 854

Download in other formats: