Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_sechiba/slowproc.f90 @ 5816

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

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

Download in other formats: