Context Navigation

source: branches/publications/ORCHILEAK-Gommet/src_sechiba/slowproc.f90 @ 7746

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

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

Download in other formats: