Context Navigation

slowproc.f90 @ 6737

Last change on this file since 6737 was 5252, checked in by chao.yue, 6 years ago
commit a flag for transition involvign bioenergy
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 259.2 KB

Line
1	! =================================================================================================================================
2	! MODULE : slowproc
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Groups the subroutines that: (1) initialize all variables used in
10	!! slowproc_main, (2) prepare the restart file for the next simulation, (3) Update the
11	!! vegetation cover if needed, and (4) handle all slow processes if the carbon
12	!! cycle is activated (call STOMATE) or update the vegetation properties (LAI and
13	!! fractional cover) in the case of a run with only SECHIBA.
14	!!
15	!!\n DESCRIPTION: None
16	!!
17	!! RECENT CHANGE(S): None
18	!!
19	!! REFERENCE(S) :
20	!!
21	!! SVN :
22	!! $HeadURL$
23	!! $Date$
24	!! $Revision$
25	!! \n
26	!_ ================================================================================================================================
27
28	MODULE slowproc
29
30	USE defprec
31	USE constantes
32	USE constantes_soil
33	USE constantes_soil_var !! crop irrigation needs, xuhui
34	USE pft_parameters
35	USE ioipsl
36	USE xios_orchidee
37	USE ioipsl_para
38	USE sechiba_io_p
39	USE interpol_help
40	USE stomate
41	USE stomate_data
42	USE grid
43	USE time, ONLY : dt_sechiba, dt_stomate, one_day, FirstTsYear, LastTsDay, LastTsYear
44	USE time, ONLY : year_start, month_start, day_start, sec_start
45	USE time, ONLY : month_end, day_end
46	USE mod_orchidee_para
47
48	IMPLICIT NONE
49
50	! Private & public routines
51
52	PRIVATE
53	PUBLIC slowproc_main, slowproc_clear, slowproc_initialize, slowproc_finalize, slowproc_change_frac, slowproc_veget
54
55	!
56	! variables used inside slowproc module : declaration and initialisation
57	!
58	REAL(r_std), SAVE :: slope_default = 0.1
59	!$OMP THREADPRIVATE(slope_default)
60	INTEGER(i_std) , SAVE :: veget_update !! update frequency in years for landuse (nb of years)
61	!$OMP THREADPRIVATE(veget_update)
62	!spitfire
63	REAL(r_std), SAVE :: m_lightn_default = 0.02
64	!$OMP THREADPRIVATE(m_lightn_default)
65	LOGICAL, SAVE :: read_popdens
66	!$OMP THREADPRIVATE(read_popdens)
67	LOGICAL, SAVE :: read_humign
68	!$OMP THREADPRIVATE(read_humign)
69	REAL(r_std), SAVE :: popdens_default = 1.
70	!$OMP THREADPRIVATE(popdens_default)
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: m_lightn
72	!$OMP THREADPRIVATE(m_lightn)
73	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: glccNetLCC !! the land-cover-change (LCC) matrix in case a gross LCC is
74	!$OMP THREADPRIVATE(glccNetLCC)
75	!! used.
76	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: glccSecondShift
77	!$OMP THREADPRIVATE(glccSecondShift)
78	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: glccPrimaryShift
79	!$OMP THREADPRIVATE(glccPrimaryShift)
80
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: harvest_matrix
82	!$OMP THREADPRIVATE(harvest_matrix)
83	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: harvest_biomass
84	!$OMP THREADPRIVATE(harvest_biomass)
85	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: bound_spa
86	!$OMP THREADPRIVATE(bound_spa)
87	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: newvegfrac !! fraction of different MTCs that is used to guide on how to
88	!! allocate the newly created MTC in gross LUC.
89	!$OMP THREADPRIVATE(newvegfrac)
90
91	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: proxy_anidens
92	!$OMP THREADPRIVATE(proxy_anidens)
93	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: popd
94	!$OMP THREADPRIVATE(popd)
95	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: humign
96	!$OMP THREADPRIVATE(humign)
97	REAL(r_std), SAVE :: m_ba_default = 0.
98	!$OMP THREADPRIVATE(m_ba_default)
99	LOGICAL, SAVE :: read_observed_ba
100	!$OMP THREADPRIVATE(read_observed_ba)
101	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: m_observed_ba
102	!$OMP THREADPRIVATE(m_observed_ba)
103
104	REAL(r_std), SAVE :: m_cf_coarse_default = 0.
105	!$OMP THREADPRIVATE(m_cf_coarse_default)
106	LOGICAL, SAVE :: read_cf_coarse
107	!$OMP THREADPRIVATE(read_cf_coarse)
108	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: m_cf_coarse
109	!$OMP THREADPRIVATE(m_cf_coarse)
110
111	REAL(r_std), SAVE :: m_cf_fine_default = 0.
112	!$OMP THREADPRIVATE(m_cf_fine_default)
113	LOGICAL, SAVE :: read_cf_fine
114	!$OMP THREADPRIVATE(read_cf_fine)
115	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: m_cf_fine
116	!$OMP THREADPRIVATE(m_cf_fine)
117
118	REAL(r_std), SAVE :: m_ratio_default = 0.
119	!$OMP THREADPRIVATE(m_ratio_default)
120	LOGICAL, SAVE :: read_ratio
121	!$OMP THREADPRIVATE(read_ratio)
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: m_ratio
123	!$OMP THREADPRIVATE(m_ratio)
124
125	REAL(r_std), SAVE :: m_ratio_flag_default = 0.
126	!$OMP THREADPRIVATE(m_ratio_flag_default)
127	LOGICAL, SAVE :: read_ratio_flag
128	!$OMP THREADPRIVATE(read_ratio_flag)
129	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: m_ratio_flag
130	!$OMP THREADPRIVATE(m_ratio_flag)
131
132	!endspit
133	!
134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: clayfraction !! Clayfraction (0-1, unitless)
135	!$OMP THREADPRIVATE(clayfraction)
136	INTEGER, SAVE :: printlev_loc !! Local printlev in slowproc module
137	!$OMP THREADPRIVATE(printlev_loc)
138	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: laimap !! LAI map when the LAI is prescribed and not calculated by STOMATE
139	!$OMP THREADPRIVATE(laimap)
140	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilclass_default
141	!$OMP THREADPRIVATE(soilclass_default)
142	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max_new !! New year fraction of vegetation type (0-1, unitless)
143	!$OMP THREADPRIVATE(veget_max_new)
144	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_new !! New year fraction of ice+lakes+cities+... (0-1, unitless)
145	!$OMP THREADPRIVATE(frac_nobio_new)
146	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio_lastyear !! last year fraction of ice+lakes+cities+... (0-1, unitless)
147	!$OMP THREADPRIVATE(frac_nobio_lastyear)
148	INTEGER(i_std), SAVE :: lcanop !! canopy levels used for LAI
149	!$OMP THREADPRIVATE(lcanop)
150	INTEGER(i_std) , SAVE :: veget_year !! year for vegetation update
151	!$OMP THREADPRIVATE(veget_year)
152	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vegetnew_firstday !! next year fraction of vegetation type (0-1, unitless)
153	!$OMP THREADPRIVATE(vegetnew_firstday)
154	!gmjc 15Feb2016 avoid grazing wet
155	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fc_grazing
156	!$OMP THREADPRIVATE(fc_grazing)
157	!end gmjc
158
159	CONTAINS
160
161	!! ================================================================================================================================
162	!! SUBROUTINE : slowproc_initialize
163	!!
164	!>\BRIEF Initialize slowproc module and call initialization of stomate module
165	!!
166	!! DESCRIPTION : Allocate module variables, read from restart file or initialize with default values
167	!! Call initialization of stomate module.
168	!!
169	!! MAIN OUTPUT VARIABLE(S) :
170	!!
171	!! REFERENCE(S) :
172	!!
173	!! FLOWCHART : None
174	!! \n
175	!_ ================================================================================================================================
176
177	SUBROUTINE slowproc_initialize (kjit, kjpij, kjpindex, date0, &
178	IndexLand, indexveg, lalo, neighbours, &
179	resolution, contfrac, soiltile, reinf_slope, &
180	t2m, &
181	deadleaf_cover, assim_param, lai, frac_age, &
182	height, veget, frac_nobio, njsc, &
183	veget_max, totfrac_nobio, qsintmax, rest_id, &
184	rest_id_stom, hist_id_stom, tot_bare_soil, &
185	hist_id_stom_IPCC, co2_flux, fco2_lu, temp_growth, &
186	soilc_total, thawed_humidity, depth_organic_soil, heat_Zimov, &
187	f_rot_sech, altmax)
188
189
190	!! 0.1 Input variables
191	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
192	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
193	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
194	REAL(r_std),INTENT (in) :: date0 !! Initial date of what ???
195	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
196	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
197	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file identifier
198	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
199	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the points on the land map
200	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in):: indexveg !! Indices of the points on the vegetation (3D map ???)
201	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
202	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! neighbouring grid points if land.
203	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
204	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid (0-1, unitless)
205	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m !! 2 m air temperature (K)
206	LOGICAL,DIMENSION(kjpindex), INTENT(in) :: f_rot_sech !! whether a grid is under rotation
207
208	!! 0.2 Output variables
209	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux per average ground area (gC m^{-2} dt_stomate^{-1})
210	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! CO2 flux from land-use (without forest management) (gC m^{-2} dt_stomate^{-1})
211	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! Growth temperature (°C) - Is equal to t2m_month
212	INTEGER(i_std), DIMENSION(kjpindex), INTENT(out) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
213	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (out) :: heat_Zimov !! heating associated with decomposition
214	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Leaf area index (m^2 m^{-2})
215	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height of vegetation (m)
216	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(out):: frac_age !! Age efficacity from STOMATE for isoprene
217	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type including none biological fraction (unitless)
218	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
219	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
220	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: altmax
221	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_bare_soil !! Total evaporating bare soil fraction
222	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
223	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile (0-1, unitless)
224	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: reinf_slope !! slope coef for reinfiltration
225	REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (out):: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
226	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless)
227	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water storage on vegetation from interception (mm)
228
229	!! 0.3 Modified variables
230	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (inout) :: soilc_total !! total soil carbon for use in thermal
231	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: thawed_humidity!! specified humidity of thawed soil
232	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: depth_organic_soil !! how deep is the organic soil?
233
234	!_ ================================================================================================================================
235
236	!! 1. Perform the allocation of all variables, define some files and some flags.
237	! Restart file read for Sechiba.
238	CALL slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
239	rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, reinf_slope, &
240	veget_max, tot_bare_soil, njsc, &
241	height, lcanop, veget_update, veget_year, f_rot_sech)
242
243
244	!! 2. Define Time step in days for stomate
245	dt_days = dt_stomate / one_day
246
247
248	!! 3. check time step coherence between slow processes and fast processes
249	IF ( dt_stomate .LT. dt_sechiba ) THEN
250	WRITE(numout,*) 'slow_processes: time step smaller than forcing time step, dt_sechiba=',dt_sechiba,' dt_stomate=',dt_stomate
251	CALL ipslerr_p(3,'slowproc_initialize','Coherence problem between dt_stomate and dt_sechiba',&
252	'Time step smaller than forcing time step','')
253	ENDIF
254
255	!! 4. Call stomate to initialize all variables manadged in stomate,
256	IF ( ok_stomate ) THEN
257
258	CALL stomate_initialize &
259	(kjit, kjpij, kjpindex, &
260	rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
261	indexLand, lalo, neighbours, resolution, &
262	contfrac, totfrac_nobio, clayfraction, t2m, &
263	lai, veget, veget_max, &
264	co2_flux, fco2_lu, &
265	deadleaf_cover, assim_param, thawed_humidity, depth_organic_soil, &
266	soilc_total, heat_Zimov, temp_growth, altmax)
267	ENDIF
268
269	!! 5. Specific run without the carbon cycle (STOMATE not called):
270	!! Need to initialize some variables that will be used in SECHIBA:
271	!! height, deadleaf_cover, assim_param, qsintmax.
272	IF (.NOT. ok_stomate ) THEN
273	CALL slowproc_derivvar (kjpindex, veget, lai, &
274	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
275	ELSE
276	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
277	qsintmax(:,1) = zero
278	ENDIF
279
280	END SUBROUTINE slowproc_initialize
281
282
283	!! ================================================================================================================================
284	!! SUBROUTINE : slowproc_main
285	!!
286	!>\BRIEF Main routine that manage variable initialisation (slowproc_init),
287	!! prepare the restart file with the slowproc variables, update the time variables
288	!! for slow processes, and possibly update the vegetation cover, before calling
289	!! STOMATE in the case of the carbon cycle activated or just update LAI (and possibly
290	!! the vegetation cover) for simulation with only SECHIBA
291	!!
292	!!
293	!! DESCRIPTION : (definitions, functional, design, flags): The subroutine manages
294	!! diverses tasks:
295	!! (1) Initializing all variables of slowproc (first call)
296	!! (2) Preparation of the restart file for the next simulation with all prognostic variables
297	!! (3) Compute and update time variable for slow processes
298	!! (4) Update the vegetation cover if there is some land use change (only every years)
299	!! (5) Call STOMATE for the runs with the carbone cycle activated (ok_stomate) and compute the respiration
300	!! and the net primary production
301	!! (6) Compute the LAI and possibly update the vegetation cover for run without STOMATE
302	!!
303	!! RECENT CHANGE(S): None
304	!!
305	!! MAIN OUTPUT VARIABLE(S): ::co2_flux, ::fco2_lu, ::lai, ::height, ::veget, ::frac_nobio,
306	!! ::veget_max, ::totfrac_nobio, ::soiltype, ::assim_param, ::deadleaf_cover, ::qsintmax,
307	!! and resp_maint, resp_hetero, resp_growth, npp that are calculated and stored
308	!! in stomate is activated.
309	!!
310	!! REFERENCE(S) : None
311	!!
312	!! FLOWCHART :
313	! \latexonly
314	! \includegraphics(scale=0.5){SlowprocMainFlow.eps} !PP to be finalize!!)
315	! \endlatexonly
316	!! \n
317	!_ ================================================================================================================================
318
319	SUBROUTINE slowproc_main (kjit, kjpij, kjpindex, date0, &
320	IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltile, &
321	t2m, temp_sol, stempdiag, &
322	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, &
323	!spitfire
324	wspeed, &
325	!endspit
326	gpp, &
327	deadleaf_cover, &
328	assim_param, &
329	lai, frac_age, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
330	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
331	co2_flux, fco2_lu, temp_growth,&
332	swdown, evapot_corr, & ! crops, xuhui
333	tdeep, hsdeep_long, snow, heat_Zimov, pb, &
334	sfluxCH4_deep, sfluxCO2_deep, &
335	thawed_humidity, depth_organic_soil, zz_deep, zz_coef_deep, &
336	soilc_total,snowdz,snowrho, &
337	tot_bare_soil, f_rot_sech, rot_cmd, &
338	!gmjc
339	tmc_topgrass, humcste_use,altmax)
340	! end gmjc
341
342	!! INTERFACE DESCRIPTION
343
344	!! 0.1 Input variables
345
346	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
347	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of the un-compressed grid
348	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
349	REAL(r_std),INTENT (in) :: date0 !! Initial date of what ???
350	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
351	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
352	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
353	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _history_ file identifier
354	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file identifier
355	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the points on the land map
356	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in):: indexveg !! Indices of the points on the vegetation (3D map ???)
357	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
358	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in) :: neighbours !! neighbouring grid points if land
359	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
360	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid (0-1, unitless)
361	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Relative humidity ("moisture stress") (0-1, unitless)
362	REAL(r_std), DIMENSION(kjpindex), INTENT(in) :: t2m !! 2 m air temperature (K)
363	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol !! Surface temperature (K)
364	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag !! Soil temperature (K)
365	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: shumdiag !! Relative soil moisture (0-1, unitless)
366	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: litterhumdiag !! Litter humidity (0-1, unitless)
367	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation (mm dt_stomate^{-1})
368	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation (mm dt_stomate^{-1})
369	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: gpp !! GPP of total ground area (gC m^{-2} time step^{-1}).
370	!! Calculated in sechiba, account for vegetation cover and
371	!! effective time step to obtain gpp_d
372	!!!!! crops
373	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !!downward shortwave radiation
374	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !!potential evaportranspiration (mm)
375	LOGICAL, DIMENSION(kjpindex), INTENT(out) :: f_rot_sech
376	INTEGER(i_std), DIMENSION(kjpindex,rot_cmd_max), INTENT(out) :: rot_cmd
377	!!!!! crops, xuhui
378	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: tdeep !! deep temperature profile
379	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (in) :: hsdeep_long!! deep long term soil humidity profile
380	REAL(r_std), DIMENSION(kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2]
381	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure
382	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_deep !! deep vertical profile
383	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_coef_deep!! deep vertical profile
384	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in) :: snowdz !! snow depth for each layer
385	REAL(r_std), DIMENSION(kjpindex,nsnow),INTENT(in) :: snowrho !! snow density for each layer
386	!spitfire
387	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: wspeed !!Wind speed (m/s)
388	!endspit
389	!gmjc top 5 layer grassland soil moisture for grazing
390	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tmc_topgrass
391	!end gmjc
392	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: humcste_use
393	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(inout) :: altmax
394	!! 0.2 Output variables
395	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: co2_flux !! CO2 flux per average ground area (gC m^{-2} dt_stomate^{-1})
396	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fco2_lu !! CO2 flux from land-use (without forest management) (gC m^{-2} dt_stomate^{-1})
397	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! Growth temperature (°C) - Is equal to t2m_month
398	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (out) :: heat_Zimov !! heating associated with decomposition
399	REAL(r_std), DIMENSION(kjpindex), INTENT (out) :: sfluxCH4_deep !! surface flux of CH4 to atmosphere from permafrost
400	REAL(r_std), DIMENSION(kjpindex), INTENT (out) :: sfluxCO2_deep !! surface flux of CO2 to atmosphere from permafrost
401	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
402
403	!! 0.3 Modified variables
404	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: lai !! Leaf area index (m^2 m^{-2})
405	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: height !! height of vegetation (m)
406	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(inout):: frac_age !! Age efficacity from STOMATE for isoprene
407	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: veget !! Fraction of vegetation type including none biological fractionin the mesh (unitless)
408	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
409	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
410	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
411	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(inout) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
412	REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (inout):: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
413	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless)
414	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintmax !! Maximum water storage on vegetation from interception (mm)
415	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: thawed_humidity !! specified humidity of thawed soil
416	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: depth_organic_soil !! how deep is the organic soil?
417	REAL(r_std), DIMENSION(kjpindex,ndeep,nvm), INTENT (inout) :: soilc_total !! total soil carbon for use in thermal
418
419	!! 0.4 Local variables
420	INTEGER(i_std) :: j, jv, ji !! indices
421	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_maint !! Maitanance component of autotrophic respiration in (gC m^{-2} dt_stomate^{-1})
422	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_hetero !! heterotrophic resp. (gC/(m**2 of total ground)/time step)
423	REAL(r_std), DIMENSION(kjpindex,nvm) :: resp_growth !! Growth component of autotrophic respiration in gC m^{-2} dt_stomate^{-1})
424	REAL(r_std), DIMENSION(kjpindex,nvm) :: npp !! Net Ecosystem Exchange (gC/(m**2 of total ground)/time step)
425	REAL(r_std),DIMENSION (kjpindex,nvm) :: veget_nextyear !! Temporary variable for new veget_max
426	REAL(r_std),DIMENSION (kjpindex,nnobio) :: frac_nobio_nextyear !! Temporary variable for new frac_nobio
427	REAL(r_std),DIMENSION (kjpindex) :: totfrac_nobio_lastyear !! Total fraction for the previous year
428	REAL(r_std),DIMENSION (kjpindex) :: totfrac_nobio_new !! Total fraction for the next year
429	!spitfire
430	REAL(r_std),DIMENSION (kjpindex) :: lightn !!lightning (flashes/km2/day)
431	REAL(r_std),DIMENSION (kjpindex) :: animal_density
432	REAL(r_std),DIMENSION (kjpindex) :: observed_ba !!observed burned area (ha)
433	INTEGER(i_std) :: yy, mm, dd
434	REAL(r_std) :: ss
435	REAL(r_std),DIMENSION (kjpindex) :: cf_coarse !!observed burned area (ha)
436	REAL(r_std),DIMENSION (kjpindex) :: cf_fine !!observed burned area (ha)
437	REAL(r_std),DIMENSION (kjpindex) :: ratio !!observed burned area (ha)
438	REAL(r_std),DIMENSION (kjpindex) :: ratio_flag !!observed burned area (ha)
439	REAL(r_std) :: cropshare_old, cropshare_new
440	!endspit
441	INTEGER(i_std) :: ivma,ivm,jvm !! Indices
442	INTEGER(i_std) , SAVE :: veget_year_tmp !! year for landuse
443
444	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: manage !! temporary matrix for rotation input
445	INTEGER(i_std) :: yrlen
446	CHARACTER(LEN=30) :: strManage, strVar
447	LOGICAL :: is_update_rotation_cycle = .FALSE. !! enable update of rotation cycle for CROP module
448
449	REAL(r_std) :: in_julian
450	!_ ================================================================================================================================
451
452	!! 1. Compute and update all variables linked to the date and time
453	IF (printlev_loc>=5) WRITE(numout,*) 'Entering slowproc_main' ! year_start, month_start, day_start, sec_start='
454	! year_start, month_start,day_start,sec_start
455
456	! IF ( sec_start == dt_sechiba .AND. month==1 .AND. day==1 ) THEN
457	! The current time step is the first sechiba timestep of a new year
458	! IF (printlev_loc>=4) WRITE(numout,*) "This is a new day and a new year: month, day, sec_start=", month, day, sec_start
459	! FirstTsYear=.TRUE.
460	! ELSE
461	! FirstTsYear=.FALSE.
462	! END IF
463
464	! IF ( sec_start == 0 ) THEN
465	! The current time step is the last sechiba time step on a day
466	! LastTsDay=.TRUE.
467
468	! IF ( month == 1 .AND. day == 1 ) THEN
469	! The current time step is the last sechiba time step on a year
470	! JG : note that month=1, day=1, sec=O is the last day of the year.
471	! This is due to a problem before the first call to slowproc.
472	! slowproc_main enters the first time on the 2nd time step (1800s)
473	! LastTsYear = .TRUE.
474	! IF (printlev_loc>=4) WRITE(numout,*) "This is the last sechiba time step of a year, LastTsYear is activated"
475	! ELSE
476	! LastTsYear = .FALSE.
477	! END IF
478	! ELSE
479	! LastTsDay = .FALSE.
480	! LastTsYear = .FALSE.
481	! END IF
482
483	!! 2. Activate slow processes if it is the end of the day
484	IF ( LastTsDay ) THEN
485	! 3.2.2 Activate slow processes in the end of the day
486	do_slow = .TRUE.
487
488	! 3.2.3 Count the number of days
489	days_since_beg = days_since_beg + 1
490	IF (printlev_loc>=4) WRITE(numout,*) "New days_since_beg : ",days_since_beg
491	ELSE
492	do_slow = .FALSE.
493	ENDIF
494
495	!! 3. Update the vegetation if it is time to do so.
496	!! This is done at the first sechiba time step on a new year and only every "veget_update" years.
497	!! veget_update correspond to a number of years between each vegetation updates.
498	!! Nothing is done if veget_update=0.
499	!! Update of the vegetation map can not be done if map_pft_format=false.
500
501	IF ( (map_pft_format) .AND. (veget_update .GT. 0) ) THEN
502
503	IF ( days_since_beg == 1 ) THEN !!! This is in case the simulation did not start from the beginning of the year (noted xuhui)
504	!
505	veget_year_tmp = veget_year + 1
506
507	! Update of the vegetation cover with Land Use only if
508	! the current year match the requested condition (a multiple of
509	! "veget_update")
510	IF ( MOD(veget_year_tmp - veget_year_orig, veget_update) == 0 ) THEN
511
512	WRITE(numout,*) 'We read the new vegetmax map for year =' , veget_year_tmp
513
514	! Call the routine to update the vegetation (output is vegetnew_firstday)
515	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
516	& veget_max, veget_nextyear, frac_nobio_nextyear, veget_year, .FALSE.)
517	!!!! veget_nextyear is saved for used next year, but since it is
518	!now input veget_map in the FirstTsYear, it is useless... (noted xuhui)
519	!print *,'veget_nextyear after reading the new map firsttime zz',veget_nextyear
520	!
521	ENDIF
522	!
523	ENDIF
524	ENDIF
525
526	! Update vegetation and fraction and save old values
527	frac_nobio_lastyear = frac_nobio
528
529	IF ( map_pft_format .AND. (veget_update > 0) .AND. FirstTsYear ) THEN
530	veget_year = veget_year + 1
531
532	! Update of the vegetation cover with Land Use only if
533	! the current year match the requested condition (a multiple of "veget_update")
534	IF ( MOD(veget_year - veget_year_orig, veget_update) == 0 ) THEN
535	IF (printlev_loc>=1) WRITE(numout,*) 'We are updating the vegetation map for year =' , veget_year
536
537	! Read the new the vegetation from file. Output is veget_nextyear and frac_nobio_nextyear.
538	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
539	veget_max, veget_nextyear, frac_nobio_nextyear, veget_year, .FALSE.)
540
541	IF (.NOT. use_age_class) THEN
542
543	IF (ok_rotate) THEN
544	! when rotation is activated, the conversion among different croplands will
545	! not follow the vegetation map but follow the rotation commands.
546	! Vegetation map is used as change of cropland share in the grid.
547	! This change is shared proportionally to contemporary croplands.
548	! xuhui 20160503
549	DO ji = 1,kjpindex
550	cropshare_old = SUM(veget_max(ji,:),MASK=ok_LAIdev(:))
551	cropshare_new = SUM(veget_nextyear(ji,:),MASK=ok_LAIdev(:))
552	IF (printlev>=4) THEN
553	WRITE(numout,*) 'ji, cropshare_old, cropshare_new',ji, cropshare_old, cropshare_new
554	ENDIF
555	IF (cropshare_old .LT. min_sechiba) THEN
556	! special case 1, no croplands previously existed
557	IF (printlev>=4) WRITE(numout,*) 'new croplands included, no change to veget_nextyear'
558	ELSEIF (cropshare_new .LT. min_sechiba) THEN
559	!special case 2, all croplands are removed
560	IF (printlev>=4) WRITE(numout,*) 'all croplands are killed, a case has not been tested'
561	!xuhui: need to consider how could it be well treated
562	ELSE
563	DO jv = 2,nvm
564	IF (ok_LAIdev(jv)) THEN
565	veget_nextyear(ji,jv) = veget_max(ji,jv) * cropshare_new / cropshare_old
566	ENDIF
567	ENDDO
568	ENDIF
569	ENDDO
570	ENDIF
571	veget_max_new = veget_nextyear
572	frac_nobio_new = frac_nobio_nextyear
573	! [chaoyue]
574	! the new veget_max will be calculated in the gross land use change
575	! module when use_age_class is True. Here we will just make the
576	! veget_max_new and frac_nobio_new being the old values.
577
578	!! !!!! xuhui: so the veget_max no longer updates here, so re-consider
579	!! !how to maintain crop fraction....
580	!! !!! should maintain veget_max_new so that when slowproc_change_frac
581	!! !applies, it works fine
582	!! IF (ok_rotate) THEN
583	!! DO jv = 1,nvm
584	!! !!!! proportionally keeping the bare soil and croplands as it
585	!! !was previously in order to maintain the rotation cycle
586	!! ENDDO
587	!! ENDIF
588	ELSE !!!! if use_age_class, readvegetmax is no longer used
589	veget_max_new = veget_max
590	frac_nobio_new = frac_nobio_lastyear
591	ENDIF
592
593	! Verification and correction on veget_max, calculation of veget and soiltile.
594	! [chaoyue] this call in the trunk equivalent is lacking, it's done
595	! in sechiba.f90 by calling slowproc_change_frac
596	CALL slowproc_veget (kjpindex, f_rot_sech, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile)
597
598	! Set the flag do_now_stomate_lcchange to activate stomate_lcchange.
599	! This flag will be kept to true until stomate_lcchange has been done.
600	! The variable totfrac_nobio_lastyear will only be used in stomate when this flag is activated
601	do_now_stomate_lcchange=.TRUE.
602	IF ( .NOT. ok_stomate ) THEN
603	! Special case if stomate is not activated : set the variable done_stomate_lcchange=true
604	! so that the subroutine slowproc_change_frac will be called in the end of sechiba_main.
605	done_stomate_lcchange=.TRUE.
606	END IF
607
608	!print *,'veget_max after reading the new map zz',veget_max
609	!print *,'frac_nobio_lastyear after reading the new map zz',frac_nobio_lastyear
610	!print *,'totfrac_nobio after reading the new map zz',totfrac_nobio
611	ENDIF
612	ENDIF
613
614	!WRITE(numout,*),'do_now_stomate_lcchange in slowproc,',do_now_stomate_lcchange
615	!WRITE(numout,*),'veget_year in slowproc,',veget_year
616	IF ( (rotation_update .GT. 0) .AND. FirstTsYear ) THEN
617	! rotation_update is necessarily zero, if not ok_rotate
618	IF ( MOD(veget_year - veget_year_orig, rotation_update) == 0 ) THEN ! update rotation cycle
619	IF (printlev_loc>=1) THEN
620	WRITE(numout,*) 'xuhui: updating rotation system at year ', veget_year
621	WRITE(numout,*) 'rotation_update, ', rotation_update
622	ENDIF
623	is_update_rotation_cycle = .TRUE. ! to stomate_main
624	ENDIF ! start rotation update
625	ENDIF ! FirstTsYear
626
627
628	! 5 call STOMATE, either because we want to keep track of
629	! long-term variables (WATCHOUT case) or just because STOMATE is
630	! activated
631
632	!spitfire
633	in_julian = itau2date(kjit, date0, dt_sechiba)
634	CALL ju2ymds(in_julian, yy, mm, dd, ss)
635	lightn(:)=m_lightn(:,mm)
636	animal_density(:)=proxy_anidens(:,mm)
637	observed_ba(:)=m_observed_ba(:,mm)
638	cf_coarse(:)=m_cf_coarse(:,mm)
639	cf_fine(:)=m_cf_fine(:,mm)
640	ratio(:)=m_ratio(:,mm)
641	ratio_flag(:)=m_ratio_flag(:,mm)
642	!endspit
643
644	IF ( ok_stomate ) THEN
645	! Caluclate totfrac_nobio_lastyear
646	totfrac_nobio_lastyear(:) = zero
647	DO jv = 1, nnobio
648	totfrac_nobio_lastyear(:) = totfrac_nobio_lastyear(:) + frac_nobio_lastyear(:,jv)
649	ENDDO
650
651	! Caluclate totfrac_nobio_new only for the case when the land use map has been read previously
652	IF (do_now_stomate_lcchange) THEN
653	IF (.NOT. use_age_class) THEN
654	totfrac_nobio_new(:) = zero
655	DO jv = 1, nnobio
656	totfrac_nobio_new(:) = totfrac_nobio_new(:) + frac_nobio_new(:,jv)
657	ENDDO
658	ELSE
659	totfrac_nobio_new = totfrac_nobio_lastyear
660	ENDIF
661	ELSE
662	totfrac_nobio_new(:) = zero
663	END IF
664
665	!! 4.1 Call stomate main routine that will call all c-cycle routines !
666	CALL stomate_main (kjit, kjpij, kjpindex, &
667	IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio_lastyear, clayfraction, &
668	t2m, temp_sol, stempdiag, &
669	humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, &
670	!spitfire
671	wspeed, lightn, popd, read_observed_ba, observed_ba, humign, &
672	read_cf_fine,cf_fine,read_cf_coarse,cf_coarse,read_ratio_flag,ratio_flag,read_ratio,ratio,&
673	!endspit
674	gpp, &
675	deadleaf_cover, &
676	assim_param, &
677	lai, frac_age, height, veget, veget_max, &
678	veget_max_new,vegetnew_firstday, totfrac_nobio_new, &
679	glccNetLCC,glccSecondShift,glccPrimaryShift, &
680	harvest_matrix, harvest_biomass,bound_spa, newvegfrac, &
681	hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
682	co2_flux, fco2_lu, resp_maint,resp_hetero,resp_growth,temp_growth,&
683	swdown, evapot_corr, is_update_rotation_cycle, & !!! xuhui added for crops
684	tdeep, hsdeep_long, snow, heat_Zimov, pb, &
685	sfluxCH4_deep, sfluxCO2_deep, &
686	thawed_humidity, depth_organic_soil, zz_deep, &
687	zz_coef_deep, soilc_total,snowdz,snowrho, &
688	LastTsYear, f_rot_sech, rot_cmd, &
689	!gmjc top 5 layer grassland soil moisture for grazing
690	tmc_topgrass,fc_grazing,humcste_use,altmax)
691	!end gmjc
692
693	! [chaoyue] in case of use_age_class and gross land use change, veget_max
694	! will be updated in the land use change module rather than using veget_max_new,
695	! thus we should pass the updated veget_max into veget_max_new
696	IF (done_stomate_lcchange .AND. use_age_class) THEN
697	veget_max_new = veget_max
698
699	IF (min_vegfrac .GT. min_stomate) THEN
700	WRITE(numout,*) "the vlaue of min_vegfrac is too big: ",min_vegfrac
701	WRITE(numout,*) "When calling slowproc_veget, a too small vlaue of"
702	WRITE(numout,*) "min_vegfrac will leads to readjustment of veget_max"
703	WRITE(numout,*) "and will break the conservation of veget_max in land use change"
704	STOP
705	ENDIF
706	ENDIF
707	!! 4.2 Output the respiration terms and the net primary
708	!! production (NPP) that are calculated in STOMATE
709
710	! 4.2.1 Output the 3 respiration terms
711	CALL xios_orchidee_send_field("maint_resp",resp_maint/dt_sechiba)
712	CALL xios_orchidee_send_field("hetero_resp",resp_hetero/dt_sechiba)
713	CALL xios_orchidee_send_field("growth_resp",resp_growth/dt_sechiba)
714
715	CALL histwrite_p(hist_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
716	CALL histwrite_p(hist_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg)
717	CALL histwrite_p(hist_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg)
718
719	! 4.2.2 Compute the net primary production as the diff from
720	! Gross primary productin and the growth and maintenance
721	! respirations
722	npp(:,1)=zero
723	DO j = 2,nvm
724	npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j)
725	ENDDO
726
727	CALL xios_orchidee_send_field("npp",npp/dt_sechiba)
728
729	CALL histwrite_p(hist_id, 'npp', kjit, npp, kjpindex*nvm, indexveg)
730
731	IF ( hist2_id > 0 ) THEN
732	CALL histwrite_p(hist2_id, 'maint_resp', kjit, resp_maint, kjpindex*nvm, indexveg)
733	CALL histwrite_p(hist2_id, 'hetero_resp', kjit, resp_hetero, kjpindex*nvm, indexveg)
734	CALL histwrite_p(hist2_id, 'growth_resp', kjit, resp_growth, kjpindex*nvm, indexveg)
735	CALL histwrite_p(hist2_id, 'npp', kjit, npp, kjpindex*nvm, indexveg)
736	ENDIF
737
738	ELSE
739	!! ok_stomate is not activated
740	!! Define the CO2 flux from the grid point to zero (no carbone cycle)
741	co2_flux(:,:) = zero
742	ENDIF
743
744
745	!! 5. Do daily processes if necessary
746	!!
747	IF ( do_slow ) THEN
748
749	!! 5.1 Calculate the LAI if STOMATE is not activated
750	IF ( .NOT. ok_stomate ) THEN
751	CALL slowproc_lai (kjpindex, lcanop,stempdiag, &
752	lalo,resolution,lai,laimap)
753
754	frac_age(:,:,1) = un
755	frac_age(:,:,2) = zero
756	frac_age(:,:,3) = zero
757	frac_age(:,:,4) = zero
758	ENDIF
759
760	! !! 5.2.0 crop rotation, if rotation started
761	! DO ji = 1,kjpindex
762	! DO jv = 2,nvm
763	! IF (ok_LAIdev(jv) .AND. f_rot_sechiba(ji,jv)) THEN !! whether to rotate
764	!
765	! ENDIF
766	! ENDDO
767	! ENDDO
768	!
769	! !! end rotation, xuhui
770
771	!! 5.2 Update veget
772	CALL slowproc_veget (kjpindex, f_rot_sech, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile)
773
774	!! 5.3 updates qsintmax and other derived variables
775	IF ( .NOT. ok_stomate ) THEN
776	CALL slowproc_derivvar (kjpindex, veget, lai, &
777	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
778	ELSE
779	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
780	qsintmax(:,1) = zero
781	ENDIF
782	END IF
783
784	!! 6. Calculate tot_bare_soil needed in hydrol, diffuco and condveg (fraction in the mesh)
785	tot_bare_soil(:) = veget_max(:,1)
786	DO jv = 2, nvm
787	DO ji =1, kjpindex
788	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
789	ENDDO
790	END DO
791
792
793	!! 7. Do some basic tests on the surface fractions updated above, only if
794	!! slowproc_veget has been done (do_slow). No change of the variables.
795	IF (do_slow) THEN
796	CALL slowproc_checkveget(kjpindex, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
797	END IF
798
799	!! 8. Write output fields
800	CALL xios_orchidee_send_field("tot_bare_soil",tot_bare_soil)
801
802	IF ( .NOT. almaoutput) THEN
803	CALL histwrite_p(hist_id, 'tot_bare_soil', kjit, tot_bare_soil, kjpindex, IndexLand)
804	END IF
805
806
807	IF (printlev_loc>=3) WRITE (numout,*) ' slowproc_main done '
808
809	END SUBROUTINE slowproc_main
810
811
812	!! ================================================================================================================================
813	!! SUBROUTINE : slowproc_finalize
814	!!
815	!>\BRIEF Write to restart file variables for slowproc module and call finalization of stomate module
816	!!
817	!! DESCRIPTION :
818	!!
819	!! MAIN OUTPUT VARIABLE(S) :
820	!!
821	!! REFERENCE(S) :
822	!!
823	!! FLOWCHART : None
824	!! \n
825	!_ ================================================================================================================================
826
827	SUBROUTINE slowproc_finalize (kjit, kjpindex, rest_id, IndexLand, &
828	njsc, lai, height, veget, &
829	frac_nobio, veget_max, reinf_slope, &
830	zz_deep, zz_coef_deep, thawed_humidity, depth_organic_soil, &
831	assim_param, frac_age, altmax)
832
833	!! 0.1 Input variables
834	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
835	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
836	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
837	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the points on the land map
838	INTEGER(i_std), DIMENSION(kjpindex), INTENT(in) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
839	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! Leaf area index (m^2 m^{-2})
840	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: height !! height of vegetation (m)
841	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type including none biological fraction (unitless)
842	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: altmax
843	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
844	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
845	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: reinf_slope !! slope coef for reinfiltration
846	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_deep !! deep vertical profile
847	REAL(r_std), DIMENSION(ndeep), INTENT (in) :: zz_coef_deep !! deep vertical profile
848	REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (in):: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
849	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT(in):: frac_age !! Age efficacity from STOMATE for isoprene
850	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: thawed_humidity!! specified humidity of thawed soil
851	REAL(r_std), DIMENSION(kjpindex), INTENT (inout) :: depth_organic_soil !! how deep is the organic soil?
852
853	!! 0.4 Local variables
854	REAL(r_std) :: tmp_day(1) !! temporary variable for I/O
855	INTEGER :: jf !! Indice
856	CHARACTER(LEN=4) :: laistring !! Temporary character string
857	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
858	!_ ================================================================================================================================
859
860	IF (printlev_loc>=3) WRITE (numout,*) 'Write restart file with SLOWPROC variables '
861
862	! 2.1 Write a series of variables controled by slowproc: day
863	! counter, vegetation fraction, max vegetation fraction, LAI
864	! variable from stomate, fraction of bare soil, soiltype
865	! fraction, clay fraction, height of vegetation, map of LAI
866
867	CALL restput_p (rest_id, 'veget', nbp_glo, nvm, 1, kjit, veget, 'scatter', nbp_glo, index_g)
868
869	CALL restput_p (rest_id, 'veget_max', nbp_glo, nvm, 1, kjit, veget_max, 'scatter', nbp_glo, index_g)
870
871	CALL restput_p (rest_id, 'lai', nbp_glo, nvm, 1, kjit, lai, 'scatter', nbp_glo, index_g)
872
873	CALL restput_p (rest_id, 'frac_nobio', nbp_glo, nnobio, 1, kjit, frac_nobio, 'scatter', nbp_glo, index_g)
874
875	CALL restput_p (rest_id, 'frac_age', nbp_glo, nvm, nleafages, kjit, frac_age, 'scatter', nbp_glo, index_g)
876
877	! Add the soil_classif as suffix for the variable name of njsc when it is stored in the restart file.
878	IF (soil_classif == 'zobler') THEN
879	var_name= 'njsc_zobler'
880	ELSE IF (soil_classif == 'usda') THEN
881	var_name= 'njsc_usda'
882	END IF
883	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, REAL(njsc, r_std), 'scatter', nbp_glo, index_g)
884
885	IF ( hydrol_cwrr ) THEN
886	CALL restput_p (rest_id, 'reinf_slope', nbp_glo, 1, 1, kjit, reinf_slope, 'scatter', nbp_glo, index_g)
887	END IF
888
889	CALL restput_p (rest_id, 'clay_frac', nbp_glo, 1, 1, kjit, clayfraction, 'scatter', nbp_glo, index_g)
890	!
891	! The height of the vegetation could in principle be recalculated at the beginning of the run.
892	! However, this is very tedious, as many special cases have to be taken into account. This variable
893	! is therefore saved in the restart file.
894	CALL restput_p (rest_id, 'height', nbp_glo, nvm, 1, kjit, height, 'scatter', nbp_glo, index_g)
895	!
896	! Specific case where the LAI is read and not calculated by STOMATE: need to be saved
897	IF (read_lai) THEN
898	CALL restput_p (rest_id, 'laimap', nbp_glo, nvm, 12, kjit, laimap)
899	ENDIF
900	!
901	! If there is some land use change, write the year for the land use ???
902	IF (map_pft_format) THEN
903	CALL restput_p (rest_id, 'veget_year', kjit, veget_year)
904	ENDIF
905
906	! 2.2 Write restart variables managed by STOMATE
907	IF ( ok_stomate ) THEN
908	CALL stomate_finalize (kjit, kjpindex, indexLand, clayfraction, &
909	zz_deep, zz_coef_deep, thawed_humidity, depth_organic_soil, &
910	assim_param,altmax)
911	ENDIF
912
913	END SUBROUTINE slowproc_finalize
914
915
916	!! ================================================================================================================================
917	!! SUBROUTINE : slowproc_init
918	!!
919	!>\BRIEF Initialisation of all variables linked to SLOWPROC
920	!!
921	!! DESCRIPTION : (definitions, functional, design, flags): The subroutine manages
922	!! diverses tasks:
923	!!
924	!! RECENT CHANGE(S): None
925	!!
926	!! MAIN OUTPUT VARIABLE(S): ::lcanop, ::veget_update, ::veget_year,
927	!! ::lai, ::veget, ::frac_nobio, ::totfrac_nobio, ::veget_max, ::height, ::soiltype
928	!!
929	!! REFERENCE(S) : None
930	!!
931	!! FLOWCHART : None
932	!! \n
933	!_ ================================================================================================================================
934
935	SUBROUTINE slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
936	rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, reinf_slope, &
937	veget_max, tot_bare_soil, njsc, &
938	height, lcanop, veget_update, veget_year, f_rot_sech)
939
940	!! INTERFACE DESCRIPTION
941
942	!! 0.1 Input variables
943	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
944	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - Terrestrial pixels only
945	INTEGER(i_std), INTENT (in) :: rest_id !! Restart file identifier
946
947	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: IndexLand !! Indices of the land points on the map
948	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
949	INTEGER(i_std), DIMENSION (kjpindex,NbNeighb), INTENT(in):: neighbours !! Vector of neighbours for each grid point
950	!! (1=North and then clockwise)
951	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x and y of the grid (m)
952	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid (unitless)
953
954	!! 0.2 Output variables
955	INTEGER(i_std), INTENT(out) :: lcanop !! Number of Canopy level used to compute LAI
956	INTEGER(i_std), INTENT(out) :: veget_update !! update frequency in timesteps (years) for landuse
957	INTEGER(i_std), INTENT(out) :: veget_year !! first year for landuse (year or index ???)
958
959	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: lai !! Leaf Area index (m^2 / m^2)
960	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget !! Fraction of vegetation type in the mesh (unitless)
961	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: frac_nobio !! Fraction of ice,lakes,cities, ... in the mesh (unitless)
962	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio !! Total fraction of ice+lakes+cities+... in the mesh (unitless)
963	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: veget_max !! Max fraction of vegetation type in the mesh (unitless)
964	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
965	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! Height of vegetation or surface in genral ??? (m)
966	REAL(r_std),DIMENSION (kjpindex,nvm,nleafages), INTENT (out):: frac_age !! Age efficacity from STOMATE for isoprene
967	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
968	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: reinf_slope !! slope coef for reinfiltration
969	INTEGER(i_std), DIMENSION(kjpindex), INTENT(out) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
970	LOGICAL, DIMENSION(kjpindex), INTENT(in) :: f_rot_sech !! whether a grid is under rotation
971
972	!! 0.3 Local variables
973	REAL(r_std) :: tmp_veget_year(1) !! temporary variable
974	REAL(r_std) :: zcanop !! ???? soil depth taken for canopy
975	INTEGER(i_std) :: vtmp(1) !! temporary variable
976	REAL(r_std), DIMENSION(nslm) :: zsoil !! soil depths at diagnostic levels
977	CHARACTER(LEN=4) :: laistring !! Temporary character string
978	INTEGER(i_std) :: l, jf !! Indices
979	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
980	INTEGER(i_std) :: ji, jv, ier,jst !! Indices
981	LOGICAL :: get_slope
982	REAL(r_std) :: frac_nobio1 !! temporary variable for frac_nobio(see above)
983	REAL(r_std), DIMENSION(kjpindex) :: tmp_real
984	REAL(r_std), DIMENSION(kjpindex,nslm) :: stempdiag2_bid !! matrix to store stempdiag_bid
985	REAL(r_std), DIMENSION (kjpindex,nscm) :: soilclass !! Fractions of each soil textural class in the grid cell (0-1, unitless)
986	CHARACTER(LEN=30), SAVE :: veget_str !! update frequency for landuse
987	!$OMP THREADPRIVATE(veget_str)
988	REAL(r_std), DIMENSION(kjpindex) :: frac_crop_tot !! Total fraction occupied by crops (0-1, unitless)
989	REAL(r_std),DIMENSION (kjpindex,nvm) :: veget_nextyear !! Temporary variable for new veget_max
990	REAL(r_std),DIMENSION (kjpindex,nnobio) :: frac_nobio_nextyear!! Temporary variable for new frac_nobio
991	LOGICAL :: found_restart !! found_restart=true if all 3 variables veget_max, veget and
992	!! frac_nobio are read from restart file
993
994	LOGICAL, SAVE :: read_veg_map_fr_restfile = .FALSE.
995	!spitfire
996	CHARACTER(LEN=80) :: data_filename
997	!endspit
998	INTEGER(i_std) :: ivma,ivm,jvm !! Indices
999	REAL(r_std),DIMENSION (nbp_glo,nvm) :: veget_max_g !! Fraction of vegetation type at global scale
1000	REAL(r_std),DIMENSION (nbp_glo,nnobio) :: frac_nobio_g !! Fraction of ice, lakes, cities etc. in the mesh (global)
1001	REAL(r_std),DIMENSION (kjpindex,nvmap) :: veget_max_map !! Fraction of vegetation for MTCs, used to hold the vegetation
1002	!! fractions for the first year of spin-up read from a MTC land cover map.
1003	REAL(r_std),DIMENSION (kjpindex,nvmap) :: veget_ny_map !! Fraction of vegetation for MTCs, it's a dummy variable
1004	!! used to initialization of veget_max_map for the first year
1005	!! run from scratch.
1006	!gmjc soil field capacity
1007	INTEGER(i_std) :: temp_njsc
1008	!end gmjc
1009	!_ ================================================================================================================================
1010
1011	! Initialize local printlev
1012	printlev_loc=get_printlev('slowproc')
1013	IF (printlev_loc>=3) WRITE (numout,*) "In slowproc_init"
1014
1015
1016	!! 1. Allocation
1017
1018	ALLOCATE (clayfraction(kjpindex),stat=ier)
1019	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable clayfraction','','')
1020	clayfraction(:)=undef_sechiba
1021
1022	! Initialisation of the fraction of the different vegetation: Start with 100% of bare soil
1023	ALLOCATE (soilclass_default(nscm),stat=ier)
1024	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable soilclass_default','','')
1025	soilclass_default(:)=undef_sechiba
1026
1027	! Allocation of the fraction of non biospheric areas
1028	ALLOCATE(frac_nobio_lastyear(kjpindex, nnobio), STAT=ier)
1029	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable frac_nobio_lastyear','','')
1030	frac_nobio_lastyear(:,:) = zero
1031
1032	ALLOCATE(veget_max_new(kjpindex, nvm), STAT=ier)
1033	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable veget_max_new','','')
1034
1035	! Allocation of the fraction of non biospheric areas
1036	ALLOCATE(frac_nobio_new(kjpindex, nnobio), STAT=ier)
1037	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable frac_nobio_new','','')
1038
1039	! Allocate laimap
1040	IF (read_lai)THEN
1041	ALLOCATE (laimap(kjpindex,nvm,12),stat=ier)
1042	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable laimap','','')
1043	ELSE
1044	ALLOCATE (laimap(1,1,1), stat=ier)
1045	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable laimap(1,1,1)','','')
1046	ENDIF
1047	!gmjc allocate fc_grazing
1048	ALLOCATE(fc_grazing(kjpindex), STAT=ier)
1049	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable fc_grazing','','')
1050	!end gmjc
1051	! Allocation of next year vegetation fraction in case of land use change
1052	ier=-1
1053	ALLOCATE(vegetnew_firstday(kjpindex, nvm), STAT=ier)
1054	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable vegetnew_firstday','','')
1055	vegetnew_firstday(:,1) = un
1056	vegetnew_firstday(:,2:nvm) = zero
1057
1058	!spitfire
1059	ier=-1
1060	ALLOCATE (m_lightn(kjpindex,12),stat=ier)
1061	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable m_lightn','','')
1062	m_lightn(:,:) = zero
1063	!
1064	ier=-1
1065	ALLOCATE (glccSecondShift(kjpindex,12),stat=ier)
1066	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable glccSecondShift','','')
1067	glccSecondShift(:,:) = zero
1068	!
1069	ier=-1
1070	ALLOCATE (glccPrimaryShift(kjpindex,12),stat=ier)
1071	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable glccPrimaryShift','','')
1072	glccPrimaryShift(:,:) = zero
1073	!
1074	ier=-1
1075	ALLOCATE (glccNetLCC(kjpindex,12),stat=ier)
1076	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable glccNetLCC','','')
1077	glccNetLCC(:,:) = zero
1078	!
1079	ier=-1
1080	ALLOCATE (harvest_matrix(kjpindex,12),stat=ier)
1081	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable harvest_matrix','','')
1082	harvest_matrix(:,:) = zero
1083	!
1084	ier=-1
1085	ALLOCATE (harvest_biomass(kjpindex,12),stat=ier)
1086	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable harvest_biomass','','')
1087	harvest_biomass(:,:) = zero
1088	!
1089	ier=-1
1090	ALLOCATE (bound_spa(kjpindex,nvm),stat=ier)
1091	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable bound_spa','','')
1092	bound_spa(:,:) = zero
1093	!
1094	ier=-1
1095	ALLOCATE (newvegfrac(kjpindex,nvmap),stat=ier)
1096	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable newvegfrac','','')
1097	newvegfrac(:,:) = zero
1098	!
1099	ier=-1
1100	ALLOCATE (proxy_anidens(kjpindex,12),stat=ier)
1101	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable proxy_anidens','','')
1102	proxy_anidens(:,:) = zero
1103	!
1104	ier=-1
1105	ALLOCATE (m_observed_ba(kjpindex,12),stat=ier)
1106	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable m_observed_ba','','')
1107	m_observed_ba(:,:) = zero
1108	!
1109	ier=-1
1110	ALLOCATE (m_cf_coarse(kjpindex,12),stat=ier)
1111	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable m_cf_coarse','','')
1112	m_cf_coarse(:,:) = zero
1113	!
1114	ier=-1
1115	ALLOCATE (m_cf_fine(kjpindex,12),stat=ier)
1116	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable m_cf_fine','','')
1117	m_cf_fine(:,:)= zero
1118	!
1119	ier=-1
1120	ALLOCATE (m_ratio(kjpindex,12),stat=ier)
1121	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable m_ratio','','')
1122	m_ratio(:,:) =zero
1123	!
1124	ier=-1
1125	ALLOCATE (m_ratio_flag(kjpindex,12),stat=ier)
1126	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable m_ratio_flag','','')
1127	m_ratio_flag(:,:) = zero
1128	!
1129	ier=-1
1130	ALLOCATE (popd(kjpindex),stat=ier)
1131	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable popd','','')
1132	popd(:) = zero
1133	!
1134	ier=-1
1135	ALLOCATE (humign(kjpindex),stat=ier)
1136	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable humign','','')
1137	humign(:) = zero
1138	!endspit
1139
1140	!! 2. Read variables from restart file
1141
1142	found_restart=.TRUE.
1143	var_name= 'veget'
1144	CALL ioconf_setatt_p('UNITS', '-')
1145	CALL ioconf_setatt_p('LONG_NAME','Vegetation fraction')
1146	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget, "gather", nbp_glo, index_g)
1147	IF ( ALL( veget(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
1148
1149	var_name= 'veget_max'
1150	CALL ioconf_setatt_p('UNITS', '-')
1151	CALL ioconf_setatt_p('LONG_NAME','Maximum vegetation fraction')
1152	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget_max, "gather", nbp_glo, index_g)
1153	IF ( ALL( veget_max(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
1154	!chaoprint
1155	!print *,'veget_max after reading restart',veget_max
1156
1157	! Get frac_nobio from the restart file
1158	var_name= 'frac_nobio'
1159	CALL ioconf_setatt_p('UNITS', '-')
1160	CALL ioconf_setatt_p('LONG_NAME','Special soil type fraction')
1161	CALL restget_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, .TRUE., frac_nobio, "gather", nbp_glo, index_g)
1162	IF ( ALL( frac_nobio(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
1163	!
1164	IF (map_pft_format .AND. .NOT. impveg) THEN
1165	IF (veget_reinit) THEN
1166	! Do not take the value read from restart file
1167	veget_year=veget_year_orig
1168	ELSE
1169	var_name= 'veget_year'
1170	CALL ioconf_setatt_p('UNITS', '-')
1171	CALL ioconf_setatt_p('LONG_NAME','Last year get in Land Use file.')
1172	CALL restget_p (rest_id, 'veget_year', kjit, .TRUE., val_exp, veget_year)
1173	IF (veget_year == val_exp) THEN
1174	! veget_year was not found in restart file
1175	veget_year=veget_year_orig
1176	ENDIF
1177	ENDIF
1178	! ENDIF
1179	! CALL bcast(veget_year)
1180
1181	!
1182	!Config Key = VEGET_UPDATE
1183	!Config Desc = Update vegetation frequency
1184	!Config If = MAP_PFT_FORMAT
1185	!Config Def = 0Y
1186	!Config Help = The veget datas will be update each this time step.
1187	!Config Units = [years]
1188	!
1189	veget_update=0
1190	WRITE(veget_str,'(a)') '0Y'
1191	CALL getin_p('VEGET_UPDATE', veget_str)
1192	l=INDEX(TRIM(veget_str),'Y')
1193	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
1194	IF (printlev_loc >= 2) WRITE(numout,*) "Update frequency for land use in years :",veget_update
1195
1196	! Coherence test
1197	IF (veget_update > 0 .AND. ok_dgvm .AND. .NOT. agriculture) THEN
1198	CALL ipslerr_p(3,'slowproc_init',&
1199	'The combination DGVM=TRUE, AGRICULTURE=FALSE and VEGET_UPDATE>0 is not possible', &
1200	'Set VEGET_UPDATE=0Y in run.def','')
1201	END IF
1202	ELSE
1203	! map_pft_format=FALSE or impveg=TRUE: there can not be any land use change, veget_update must be =0
1204	! Read VEGET_UPDATE from run.def and exit if it is different from 0Y
1205	veget_update=0
1206	WRITE(veget_str,'(a)') '0Y'
1207	CALL getin_p('VEGET_UPDATE', veget_str)
1208	l=INDEX(TRIM(veget_str),'Y')
1209	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
1210	IF (veget_update /= 0) THEN
1211	WRITE(numout,*) 'veget_update=',veget_update,' is not coeherent with map_pft_format=',map_pft_format,' or impveg=',impveg
1212	CALL ipslerr_p(3,'slowproc_init','Incoherent values between impveg, map_pft_format and veget_update', &
1213	'veget_update must be equal to 0 if map_pft_format=false or if impveg=true','')
1214	END IF
1215
1216	ENDIF
1217
1218	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_init : End of Land Use configuration'
1219
1220	IF ( hydrol_cwrr ) THEN
1221	var_name= 'reinf_slope'
1222	CALL ioconf_setatt_p('UNITS', '-')
1223	CALL ioconf_setatt_p('LONG_NAME','Slope coef for reinfiltration')
1224	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinf_slope, "gather", nbp_glo, index_g)
1225	END IF
1226
1227	! Below we define the soil texture of the grid-cells
1228	! Add the soil_classif as suffix for the variable name of njsc when it is stored in the restart file.
1229	IF (soil_classif == 'zobler') THEN
1230	var_name= 'njsc_zobler'
1231	ELSE IF (soil_classif == 'usda') THEN
1232	var_name= 'njsc_usda'
1233	ELSE
1234	CALL ipslerr_p(3,'slowproc_init','Non supported soil type classification','','')
1235	END IF
1236
1237	CALL ioconf_setatt_p('UNITS', '-')
1238	CALL ioconf_setatt_p('LONG_NAME','Index of soil type')
1239	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tmp_real, "gather", nbp_glo, index_g)
1240	IF ( ALL( tmp_real(:) .EQ. val_exp) ) THEN
1241	njsc (:) = undef_int
1242	ELSE
1243	njsc = NINT(tmp_real)
1244	ENDIF
1245
1246	var_name= 'clay_frac'
1247	CALL ioconf_setatt_p('UNITS', '-')
1248	CALL ioconf_setatt_p('LONG_NAME','Fraction of clay in each mesh')
1249	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., clayfraction, "gather", nbp_glo, index_g)
1250
1251	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_init : End CWRR configuration'
1252	!
1253	var_name= 'lai'
1254	CALL ioconf_setatt_p('UNITS', '-')
1255	CALL ioconf_setatt_p('LONG_NAME','Leaf area index')
1256	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., lai, "gather", nbp_glo, index_g)
1257
1258	! The height of the vegetation could in principle be recalculated at the beginning of the run.
1259	! However, this is very tedious, as many special cases have to be taken into account. This variable
1260	! is therefore saved in the restart file.
1261	var_name= 'height'
1262	CALL ioconf_setatt_p('UNITS', 'm')
1263	CALL ioconf_setatt_p('LONG_NAME','Height of vegetation')
1264	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., height, "gather", nbp_glo, index_g)
1265
1266	IF (read_lai)THEN
1267	var_name= 'laimap'
1268	CALL ioconf_setatt_p('UNITS', '-')
1269	CALL ioconf_setatt_p('LONG_NAME','Leaf area index read')
1270	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, .TRUE., laimap)
1271	ENDIF
1272
1273	CALL ioconf_setatt_p('UNITS', '-')
1274	CALL ioconf_setatt_p('LONG_NAME','Fraction of leaves in leaf age class ')
1275	CALL restget_p (rest_id, 'frac_age', nbp_glo, nvm, nleafages, kjit, .TRUE.,frac_age, "gather", nbp_glo, index_g)
1276
1277	!! 3. Some other initializations
1278
1279	!Config Key = SECHIBA_ZCANOP
1280	!Config Desc = Soil level used for canopy development (if STOMATE disactivated)
1281	!Config If = OK_SECHIBA and .NOT. OK_STOMATE
1282	!Config Def = 0.5
1283	!Config Help = The temperature at this soil depth is used to determine the LAI when
1284	!Config STOMATE is not activated.
1285	!Config Units = [m]
1286	zcanop = 0.5_r_std
1287	CALL setvar_p (zcanop, val_exp, 'SECHIBA_ZCANOP', 0.5_r_std)
1288
1289	! depth at center of the levels
1290	zsoil(1) = diaglev(1) / 2.
1291	DO l = 2, nslm
1292	zsoil(l) = ( diaglev(l) + diaglev(l-1) ) / 2.
1293	ENDDO
1294
1295	! index of this level
1296	vtmp = MINLOC ( ABS ( zcanop - zsoil(:) ) )
1297	lcanop = vtmp(1)
1298
1299	!
1300	! Interception reservoir coefficient
1301	!
1302	!Config Key = SECHIBA_QSINT
1303	!Config Desc = Interception reservoir coefficient
1304	!Config If = OK_SECHIBA
1305	!Config Def = 0.1
1306	!Config Help = Transforms leaf area index into size of interception reservoir
1307	!Config for slowproc_derivvar or stomate
1308	!Config Units = [m]
1309	CALL getin_p('SECHIBA_QSINT', qsintcst)
1310	IF (printlev >= 2) WRITE(numout, *)' SECHIBA_QSINT, qsintcst = ', qsintcst
1311
1312
1313
1314
1315	!! 4. Initialization of variables not found in restart file
1316
1317	IF ( impveg ) THEN
1318
1319	!! 4.1.a Case impveg=true: Initialization of variables by reading run.def
1320	!! The routine setvar_p will only initialize the variable if it was not found in restart file.
1321	!! We are on a point and thus we can read the information from the run.def
1322
1323	!Config Key = SECHIBA_VEGMAX
1324	!Config Desc = Maximum vegetation distribution within the mesh (0-dim mode)
1325	!Config If = IMPOSE_VEG
1326	!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
1327	!Config Help = The fraction of vegetation is read from the restart file. If
1328	!Config it is not found there we will use the values provided here.
1329	!Config Units = [-]
1330	CALL setvar_p (veget_max, val_exp, 'SECHIBA_VEGMAX', veget_ori_fixed_test_1)
1331
1332	!Config Key = SECHIBA_FRAC_NOBIO
1333	!Config Desc = Fraction of other surface types within the mesh (0-dim mode)
1334	!Config If = IMPOSE_VEG
1335	!Config Def = 0.0
1336	!Config Help = The fraction of ice, lakes, etc. is read from the restart file. If
1337	!Config it is not found there we will use the values provided here.
1338	!Config For the moment, there is only ice.
1339	!Config Units = [-]
1340	frac_nobio1 = frac_nobio(1,1)
1341	CALL setvar_p (frac_nobio1, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1)
1342	frac_nobio(:,:) = frac_nobio1
1343
1344	!Config Key = SECHIBA_LAI
1345	!Config Desc = LAI for all vegetation types (0-dim mode)
1346	!Config Def = 0., 8., 8., 4., 4.5, 4.5, 4., 4.5, 4., 2., 2., 2., 2.
1347	!Config If = IMPOSE_VEG and .NOT. STOMATE
1348	!Config Help = The maximum LAI used in the 0dim mode. The values should be found
1349	!Config in the restart file. The new values of LAI will be computed anyway
1350	!Config at the end of the current day. The need for this variable is caused
1351	!Config by the fact that the model may stop during a day and thus we have not
1352	!Config yet been through the routines which compute the new surface conditions.
1353	!Config Units = [-]
1354	! MICT: lai cannot be modified when ok_stomate is enabled. This breaks consistency between
1355	! other variables through biomass
1356	IF (.NOT. ok_stomate) CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax)
1357
1358	IF (.NOT. found_restart) THEN
1359	! Call slowproc_veget to correct veget_max and to calculate veget and soiltiles
1360	CALL slowproc_veget (kjpindex, f_rot_sech, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile)
1361	END IF
1362
1363	IF (impsoilt) THEN
1364
1365	! If njsc is not in restart file, then initialize soilclass from values
1366	! from run.def file and recalculate njsc
1367	IF ( ALL(njsc(:) .EQ. undef_int )) THEN
1368	!Config Key = SOIL_FRACTIONS
1369	!Config Desc = Fraction of the 3 soil types (0-dim mode)
1370	!Config Def = undef_sechiba
1371	!Config If = IMPOSE_VEG and IMPOSE_SOILT
1372	!Config Help = Determines the fraction for the 3 soil types
1373	!Config in the mesh in the following order : sand loam and clay.
1374	!Config Units = [-]
1375
1376	soilclass(1,:) = soilclass_default(:)
1377	CALL getin_p('SOIL_FRACTIONS',soilclass(1,:))
1378	! Assign for each grid-cell the % of the different textural classes (up to 12 if 'usda')
1379	DO ji=2,kjpindex
1380	! here we read, for the prescribed grid-cell, the % occupied by each of the soil texture classes
1381	soilclass(ji,:) = soilclass(1,:)
1382	ENDDO
1383
1384	! Simplify an heterogeneous grid-cell into an homogeneous one with the dominant texture
1385	njsc(:) = 0
1386	DO ji = 1, kjpindex
1387	! here we reduce to the dominant texture class
1388	njsc(ji) = MAXLOC(soilclass(ji,:),1)
1389	ENDDO
1390	END IF
1391
1392	!Config Key = CLAY_FRACTION
1393	!Config Desc = Fraction of the clay fraction (0-dim mode)
1394	!Config Def = 0.2
1395	!Config If = IMPOSE_VEG and IMPOSE_SOIL
1396	!Config Help = Determines the fraction of clay in the grid box.
1397	!Config Units = [-]
1398
1399	! If clayfraction was not in restart file it will be read fro run.def file instead of deduced
1400	! based on fractions of each textural class
1401	CALL setvar_p (clayfraction, val_exp, 'CLAY_FRACTION', clayfraction_default)
1402	ELSE
1403	! Case impsoilt=false and impveg=true
1404	IF ( MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp .OR. &
1405	MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int ) THEN
1406
1407	CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
1408	njsc(:) = 0
1409	DO ji = 1, kjpindex
1410	njsc(ji) = MAXLOC(soilclass(ji,:),1)
1411	ENDDO
1412	ENDIF
1413	ENDIF
1414
1415	!Config Key = REINF_SLOPE
1416	!Config Desc = Slope coef for reinfiltration
1417	!Config Def = 0.1
1418	!Config If = IMPOSE_VEG
1419	!Config Help = Determines the reinfiltration ratio in the grid box due to flat areas
1420	!Config Units = [-]
1421	!
1422	slope_default=0.0
1423	CALL setvar_p (reinf_slope, val_exp, 'SLOPE', slope_default)
1424
1425	!Config Key = SLOWPROC_HEIGHT
1426	!Config Desc = Height for all vegetation types
1427	!Config Def = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0
1428	!Config If = OK_SECHIBA
1429	!Config Help = The height used in the 0dim mode. The values should be found
1430	!Config in the restart file. The new values of height will be computed anyway
1431	!Config at the end of the current day. The need for this variable is caused
1432	!Config by the fact that the model may stop during a day and thus we have not
1433	!Config yet been through the routines which compute the new surface conditions.
1434	!Config Units = [m]
1435	CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc)
1436
1437
1438	ELSE IF ( .NOT. found_restart ) THEN
1439	!!! rotation is not concerned if no restarting of veget_max
1440
1441	!! 4.1.b Case impveg=false and no restart files: Initialization by reading vegetation map
1442
1443	! Initialize veget_max and frac_nobio
1444	IF ( map_pft_format ) THEN
1445	IF(use_age_class) THEN
1446	veget_ny_map(:,:)=zero
1447	! Note that veget_ny_map is only a dummy variable and not used at all,
1448	! because veget_ny_map takes the position of veget_last, which is used
1449	! only to retain the agriculture fraction when land cover is partially
1450	! updated.
1451	! `veget_max_map` is the variable that holds for us the MTC fraction
1452	! information from the land cover map.
1453	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
1454	& veget_ny_map, veget_max_map, frac_nobio_nextyear, veget_year, .TRUE.)
1455	! Now we need to change this into the map with age classes.
1456	! Since this is an initialization, we will just assign the whole PFT to the youngest age class.
1457	veget_max(:,:)=zero
1458	DO ivma=1,nvmap
1459	veget_max(:,start_index(ivma))=veget_max_map(:,ivma)
1460	ENDDO
1461	frac_nobio = frac_nobio_nextyear
1462
1463	!we're not using age groups so it's just "normal PFT map initialization"
1464	ELSE
1465	! Case without restart file and map_pft_format=true
1466	IF (printlev_loc>=3) WRITE(numout,*) 'Before call slowproc_readvegetmax in initialization phase without restart files'
1467	IF (printlev_loc>=3) WRITE(numout,*) 'veget_year=', veget_year
1468
1469	! Call the routine to update the vegetation (output is veget_nextyear)
1470	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
1471	veget_max, veget_nextyear, frac_nobio_nextyear, veget_year, .TRUE.)
1472	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_readvegetmax in initialization phase'
1473
1474	! Update vegetation with values read from the file
1475	veget_max = veget_nextyear
1476	frac_nobio = frac_nobio_nextyear
1477	END IF
1478
1479	ELSE
1480	! The interpolation of vegetation has changed.
1481	CALL getin_p('read_veg_map_fr_restfile', read_veg_map_fr_restfile)
1482	IF ( .NOT. read_veg_map_fr_restfile ) THEN
1483	! slowproc interpol :
1484	CALL slowproc_interpol_g(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, veget_max_g, frac_nobio_g)
1485	ELSE
1486	CALL slowproc_read_veg_restfile(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
1487	veget_max_g, frac_nobio_g)
1488	ENDIF
1489
1490	CALL scatter(veget_max_g, veget_max)
1491	CALL scatter(frac_nobio_g, frac_nobio)
1492
1493	! map_pft_format=FALSE: Read and interpolate Olson type map
1494	CALL slowproc_interpol(kjpindex, lalo, neighbours, resolution, contfrac, veget_max, frac_nobio)
1495	END IF
1496
1497	!! Reset totaly or partialy veget_max if using DGVM
1498	IF ( ok_dgvm ) THEN
1499	! If we are dealing with dynamic vegetation then all natural PFTs should be set to veget_max = 0
1500	! In case no agriculture is desired, agriculture PFTS should be set to 0 as well
1501	IF (agriculture) THEN
1502	DO jv = 2, nvm
1503	! dgvmjc consider pasture as not natural
1504	IF (natural(jv) .AND. .NOT. pasture(jv)) THEN
1505	veget_max(:,jv)=zero
1506	ENDIF
1507	ENDDO
1508
1509	! Calculate the fraction of crop for each point.
1510	! Sum only on the indexes corresponding to the non_natural pfts
1511	frac_crop_tot(:) = zero
1512	DO jv = 2, nvm
1513	! dgvmjc consider pasture as not natural
1514	IF(.NOT. natural(jv) .OR. pasture(jv)) THEN
1515	DO ji = 1, kjpindex
1516	frac_crop_tot(ji) = frac_crop_tot(ji) + veget_max(ji,jv)
1517	ENDDO
1518	ENDIF
1519	ENDDO
1520
1521	! Calculate the fraction of bare soil
1522	DO ji = 1, kjpindex
1523	veget_max(ji,1) = un - frac_crop_tot(ji) - SUM(frac_nobio(ji,:))
1524	ENDDO
1525	ELSE
1526	veget_max(:,:) = zero
1527	DO ji = 1, kjpindex
1528	veget_max(ji,1) = un - SUM(frac_nobio(ji,:))
1529	ENDDO
1530	ENDIF ! end agriculture
1531	ENDIF ! end ok_dgvm
1532
1533
1534	! correct veget_max and to calculate veget and soiltiles
1535	CALL slowproc_veget (kjpindex, f_rot_sech, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile)
1536
1537	END IF ! end impveg
1538
1539	!! 4.2 Continue initializing variables not found in restart file. Case for both impveg=true and false.
1540
1541	! Initialize laimap for the case read_lai if not found in restart file
1542	IF (read_lai) THEN
1543	IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN
1544	! Interpolation of LAI
1545	CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap)
1546	ENDIF
1547	ENDIF
1548
1549	! Initialize lai if not found in restart file and not already initialized using impveg
1550	IF ( MINVAL(lai) .EQ. MAXVAL(lai) .AND. MAXVAL(lai) .EQ. val_exp) THEN
1551	IF (read_lai) THEN
1552	stempdiag2_bid(1:kjpindex,1:nslm) = stempdiag_bid
1553	CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, &
1554	lalo,resolution,lai,laimap)
1555	ELSE
1556	! If we start from scratch, we set lai to zero for consistency with stomate
1557	lai(:,:) = zero
1558	ENDIF
1559
1560	frac_age(:,:,1) = un
1561	frac_age(:,:,2) = zero
1562	frac_age(:,:,3) = zero
1563	frac_age(:,:,4) = zero
1564	ENDIF
1565
1566	! Initialize heigth if not found in restart file and not already initialized using impveg
1567	IF ( MINVAL(height) .EQ. MAXVAL(height) .AND. MAXVAL(height) .EQ. val_exp) THEN
1568	! Impose height
1569	DO jv = 1, nvm
1570	height(:,jv) = height_presc(jv)
1571	ENDDO
1572	ENDIF
1573
1574	! Initialize clayfraction and njsc if not found in restart file and not already initialized using impveg
1575	IF ( MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp .OR. &
1576	MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int ) THEN
1577
1578	IF (printlev_loc>=4) WRITE (numout,*) 'clayfraction or njcs were not in restart file, call slowproc_soilt'
1579	CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
1580	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_soilt'
1581
1582	njsc(:) = 0
1583	DO ji = 1, kjpindex
1584	njsc(ji) = MAXLOC(soilclass(ji,:),1)
1585	ENDDO
1586	ENDIF
1587
1588	!gmjc 15Feb2016 avoid grazing wet
1589	! NOTE only for usda 12 classes soil map
1590	WRITE (numout,*) 'get fc_grazing gmjc'
1591	SELECTCASE(soil_classif)
1592	CASE('none')
1593	DO ji = 1, kjpindex
1594	temp_njsc = njsc(ji)
1595	fc_grazing(ji) = mcs_fao(temp_njsc)
1596	ENDDO
1597	CASE('zobler')
1598	DO ji = 1, kjpindex
1599	temp_njsc = njsc(ji)
1600	fc_grazing(ji) = mcs_fao(temp_njsc)
1601	ENDDO
1602	CASE("usda")
1603	DO ji = 1, kjpindex
1604	temp_njsc = njsc(ji)
1605	fc_grazing(ji) = mcs_usda(temp_njsc)
1606	! fc_grazing(ji) = 1.0
1607	ENDDO
1608	CASE DEFAULT
1609	WRITE(,) 'A non supported soil type classification has been chosen'
1610	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
1611	ENDSELECT
1612	!end gmjc
1613	!Config Key = GET_SLOPE
1614	!Config Desc = Read slopes from file and do the interpolation
1615	!Config Def = n
1616	!Config If =
1617	!Config Help = Needed for reading the slopesfile and doing the interpolation. This will be
1618	! used by the re-infiltration parametrization
1619	!Config Units = [FLAG]
1620	get_slope = .FALSE.
1621	CALL getin_p('GET_SLOPE',get_slope)
1622
1623	IF ( hydrol_cwrr ) THEN
1624	IF ( MINVAL(reinf_slope) .EQ. MAXVAL(reinf_slope) .AND. MAXVAL(reinf_slope) .EQ. val_exp .OR. get_slope) THEN
1625	IF (printlev_loc>=4) WRITE (numout,*) 'reinf_slope was not in restart file. Now call slowproc_slope'
1626
1627	CALL slowproc_slope(kjpindex, lalo, neighbours, resolution, contfrac, reinf_slope)
1628	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_slope'
1629
1630	ENDIF
1631	END IF
1632
1633
1634	!! 5. Some calculations always done, with and without restart files
1635
1636	! The variables veget, veget_max and frac_nobio were all read from restart file or initialized above.
1637	! Calculate now totfrac_nobio and soiltiles using these variables.
1638
1639	! Calculate totfrac_nobio
1640	totfrac_nobio(:) = zero
1641	DO jv = 1, nnobio
1642	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
1643	ENDDO
1644
1645	! Calculate soiltile. This variable do not need to be in the restart file.
1646	! The sum of all soiltiles makes one, and corresponds to the bio fraction
1647	! of the grid cell (called vegtot in hydrol)
1648	soiltile(:,:) = zero
1649	DO jv = 1, nvm
1650	jst = pref_soil_veg(jv)
1651	DO ji = 1, kjpindex
1652	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
1653	ENDDO
1654	ENDDO
1655	DO ji = 1, kjpindex
1656	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
1657	soiltile(ji,:)=soiltile(ji,:)/(1-totfrac_nobio(ji))
1658	ENDIF
1659	ENDDO
1660
1661	! Always calculate tot_bare_soil
1662	! Fraction of bare soil in the mesh (bio+nobio)
1663	tot_bare_soil(:) = veget_max(:,1)
1664	DO jv = 2, nvm
1665	DO ji =1, kjpindex
1666	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
1667	ENDDO
1668	END DO
1669
1670	!! 6. Verify consistency between different fractions. No change of the variables.
1671	IF (ok_stomate .AND. .NOT.disable_fire) THEN
1672	!spitfire
1673	!
1674	!Config Key = LIGHTNING
1675	!Config Desc = Read the ligntning map
1676	!Config Def = ?
1677	!Config If = NOT FIRE_DISABLE
1678	!Config Help = reads a 12 month lightning map which will
1679	!Config then be interpolated to daily values as needed.
1680	!
1681	data_filename='lightn.nc'
1682	CALL setvar_p (m_lightn, val_exp, 'm_lightn', m_lightn_default)
1683	CALL getin_p('LIGHTNING_FILE',data_filename)
1684	CALL slowproc_read_data(kjpindex, lalo, resolution, m_lightn,data_filename,'lightn')
1685
1686	!
1687	!Config Key = LCC_MATRIX
1688	!Config Desc = Read the ligntning map
1689	!Config Def = ?
1690	!Config Help = reads a 12 month lightning map which will
1691	!Config then be interpolated to daily values as needed.
1692	!
1693
1694	CALL setvar_p (glccNetLCC, val_exp, 'glccNetLCC', 0.)
1695	CALL setvar_p (glccSecondShift, val_exp, 'glccSecondShift', 0.)
1696	CALL setvar_p (glccPrimaryShift, val_exp, 'glccPrimaryShift', 0.)
1697	CALL setvar_p (harvest_matrix, val_exp, 'harvest_matrix', 0.)
1698	CALL setvar_p (harvest_biomass, val_exp, 'harvest_biomass', 0.)
1699	IF ( (use_age_class) .AND. (veget_update .GT. 0) ) THEN
1700	data_filename = 'GLUC_NET_LCC_FILE.nc'
1701	CALL getin_p('GLUC_NET_LCC_FILE',data_filename)
1702	CALL slowproc_read_data(kjpindex, lalo, resolution, glccNetLCC, data_filename, 'matrix')
1703
1704	data_filename = 'GLUC_SHIFT_SEC_FILE.nc'
1705	CALL getin_p('GLUC_SHIFT_SEC_FILE',data_filename)
1706	CALL slowproc_read_data(kjpindex, lalo, resolution, glccSecondShift, data_filename, 'matrix')
1707
1708	data_filename = 'GLUC_SHIFT_PRI_FILE.nc'
1709	CALL getin_p('GLUC_SHIFT_PRI_FILE',data_filename)
1710	CALL slowproc_read_data(kjpindex, lalo, resolution, glccPrimaryShift, data_filename, 'matrix')
1711
1712	data_filename = 'GLUC_NewVegFrac_File.nc'
1713	IF (gluc_newfrac_guide) THEN
1714	CALL getin_p('GLUC_NewVegFrac_File',data_filename)
1715	CALL slowproc_read_data(kjpindex, lalo, resolution, newvegfrac, data_filename, 'value')
1716	ENDIF
1717
1718	IF (allow_forestry_harvest) THEN
1719	data_filename = 'GLUC_FORESTRY_HARVEST_FILE.nc'
1720	CALL getin_p('GLUC_FORESTRY_HARVEST_FILE',data_filename)
1721	CALL slowproc_read_data(kjpindex, lalo, resolution, harvest_matrix, data_filename, 'matrix')
1722
1723	data_filename = 'GLUC_HARVEST_BIOMASS_FILE.nc'
1724	IF (gluc_use_harvest_biomass) THEN
1725	CALL getin_p('GLUC_HARVEST_BIOMASS_FILE',data_filename)
1726	! harvest_matrix should have the 1st dim as industrial wood havest, 2nd as fuel wood, and 3rd
1727	! dim as the fuel wood fraction. The unit should be GgC.
1728	CALL slowproc_read_data(kjpindex, lalo, resolution, harvest_biomass, data_filename, 'matrix')
1729	harvest_biomass(:,1:2) = harvest_biomass(:,1:2) * 1e9 !change to gC
1730	ENDIF
1731	ENDIF
1732
1733	IF (gluc_allow_trans_bioe) THEN
1734	data_filename = 'GLUC_TRANS_BIOE.nc'
1735	!CALL getin_p('GLUC_TRANS_BIOE1_FILE',data_filename)
1736	!CALL slowproc_read_data(kjpindex, lalo, resolution, trans_bioe1_matrix, data_filename, 'matrix')
1737	ENDIF
1738
1739	ENDIF
1740
1741	IF ( use_age_class .AND. use_bound_spa ) THEN
1742	CALL getin_p('GLUC_AGE_THRESHOLD_FILE',data_filename)
1743	CALL slowproc_read_data(kjpindex, lalo, resolution, bound_spa, data_filename, 'value')
1744	ENDIF
1745
1746	!Config Key = proxy_anidens
1747	!Config Desc = Read the general animal_density map
1748	!Config Def = ?
1749	!Config Help = reads a 12 month lightning map which will
1750	!Config then be interpolated to daily values as needed.
1751	!
1752
1753
1754	!Config Key = OBSERVED_BA_FLAG
1755	!Config Desc = Read the observed burned_area flag
1756	!Config Def = ?
1757	!
1758	read_observed_ba = .FALSE.
1759	CALL getin_p('READ_OBSERVED_BA',read_observed_ba)
1760	WRITE(numout,*) 'flag for READ_OBSERVED_BA ', read_observed_ba
1761
1762	!
1763	!Config Key = OBSERVED_BURNED_AREA
1764	!Config Desc = Read the population density map
1765	!Config Def = ?
1766	!Config Help = reads a one-year monthly burned area map with 12 as the value
1767	!of time dimension
1768	!
1769	IF(read_observed_ba) THEN
1770	CALL setvar_p (m_observed_ba, val_exp, 'm_observed_ba', m_ba_default)
1771	CALL getin_p('BA_FILE',data_filename)
1772	CALL slowproc_read_data(kjpindex, lalo, resolution, m_observed_ba,data_filename,'ba')
1773	ENDIF
1774
1775	! 1.
1776	!Config Key = OBSERVED_BA_FLAG
1777	!Config Desc = Read the observed burned_area flag
1778	!Config Def = ?
1779	!
1780	read_cf_coarse = .FALSE.
1781	CALL getin_p('READ_CF_COARSE',read_cf_coarse)
1782	WRITE(numout,*) 'flag for READ_CF_COARSE ', read_cf_coarse
1783
1784	!
1785	!Config Key = OBSERVED_BURNED_AREA
1786	!Config Desc = Read the population density map
1787	!Config Def = ?
1788	!Config Help = reads a one-year monthly burned area map with 12 as the value
1789	!of time dimension
1790	!
1791	IF(read_cf_coarse) THEN
1792	CALL setvar_p (m_cf_coarse, val_exp, 'm_cf_coarse', m_cf_coarse_default)
1793	CALL getin_p('CF_COARSE_FILE',data_filename)
1794	CALL slowproc_read_data(kjpindex, lalo, resolution, m_cf_coarse,data_filename,'cf')
1795	ENDIF
1796
1797	! 2.
1798	!Config Key = OBSERVED_BA_FLAG
1799	!Config Desc = Read the observed burned_area flag
1800	!Config Def = ?
1801	!
1802	read_cf_fine = .FALSE.
1803	CALL getin_p('READ_CF_FINE',read_cf_fine)
1804	WRITE(numout,*) 'flag for READ_CF_FINE ', read_cf_fine
1805
1806	!
1807	!Config Key = OBSERVED_BURNED_AREA
1808	!Config Desc = Read the population density map
1809	!Config Def = ?
1810	!Config Help = reads a one-year monthly burned area map with 12 as the value
1811	!of time dimension
1812	!
1813	IF(read_cf_fine) THEN
1814	CALL setvar_p (m_cf_fine, val_exp, 'm_cf_fine', m_cf_fine_default)
1815	CALL getin_p('CF_FINE_FILE',data_filename)
1816	CALL slowproc_read_data(kjpindex, lalo, resolution, m_cf_fine,data_filename,'cf')
1817	ENDIF
1818
1819	! 3.
1820	!Config Key = OBSERVED_BA_FLAG
1821	!Config Desc = Read the observed burned_area flag
1822	!Config Def = ?
1823	!
1824	read_ratio = .TRUE.
1825	CALL getin_p('READ_RATIO',read_ratio)
1826	WRITE(numout,*) 'flag for READ_RATIO ', read_ratio
1827
1828	!
1829	!Config Key = OBSERVED_BURNED_AREA
1830	!Config Desc = Read the population density map
1831	!Config Def = ?
1832	!Config Help = reads a one-year monthly burned area map with 12 as the value
1833	!of time dimension
1834	!
1835	IF(read_ratio) THEN
1836	CALL setvar_p (m_ratio, val_exp, 'm_ratio', m_ratio_default)
1837	CALL getin_p('RATIO_FILE',data_filename)
1838	CALL slowproc_read_data(kjpindex, lalo, resolution, m_ratio,data_filename,'ratio')
1839	ENDIF
1840
1841	! 4.
1842	!Config Key = OBSERVED_BA_FLAG
1843	!Config Desc = Read the observed burned_area flag
1844	!Config Def = ?
1845	!
1846	read_ratio_flag = .TRUE.
1847	CALL getin_p('READ_RATIO_FLAG',read_ratio_flag)
1848	WRITE(numout,*) 'flag for READ_RATIO_FLAG ', read_ratio_flag
1849
1850	!
1851	!Config Key = OBSERVED_BURNED_AREA
1852	!Config Desc = Read the population density map
1853	!Config Def = ?
1854	!Config Help = reads a one-year monthly burned area map with 12 as the value
1855	!of time dimension
1856	!
1857	IF(read_ratio_flag) THEN
1858	CALL setvar_p (m_ratio_flag, val_exp, 'm_ratio_flag', m_ratio_flag_default)
1859	CALL getin_p('RATIO_FLAG_FILE',data_filename)
1860	CALL slowproc_read_data(kjpindex, lalo, resolution, m_ratio_flag,data_filename,'invalid_flag')
1861	ENDIF
1862
1863	!Config Key = POPDENS_FLAG
1864	!Config Desc = Read the popdens flag
1865	!Config Def = n
1866	!
1867	read_popdens = .FALSE.
1868	CALL getin_p('READ_POPDENS',read_popdens)
1869	WRITE(numout,*) 'flag for READ_POPDENS ',read_popdens
1870
1871	!
1872	!Config Key = POPDENS
1873	!Config Desc = Read the population density map
1874	!Config Def = ?
1875	!Config Help = reads a yearly map
1876	!
1877	popd=0.
1878	IF(read_popdens) THEN
1879	CALL setvar_p (popd, val_exp, 'popdens', popdens_default)
1880	CALL getin_p('POPDENS_FILE',data_filename)
1881	CALL slowproc_read_annual(kjpindex, lalo, resolution, popd,data_filename,'popdens')
1882	ENDIF
1883
1884	!Config Key = HUMIGN_FLAG
1885	!Config Desc = Read the human ignition parameter flag
1886	!Config Def = n
1887	!
1888	read_humign = .FALSE.
1889	CALL getin_p('READ_HUMIGN',read_humign)
1890	WRITE(numout,*) 'flag for READ_HUMIGN:',read_humign
1891
1892	!
1893	!Config Key = POPDENS
1894	!Config Desc = Read the population density map
1895	!Config Def = ?
1896	!Config Help = reads a yearly map
1897	!
1898	humign=0.22
1899	IF(read_humign) THEN
1900	CALL setvar_p (humign, val_exp, 'HUMIGN_FILE', 0.22)
1901	CALL getin_p('HUMIGN_FILE',data_filename)
1902	CALL slowproc_read_annual(kjpindex, lalo, resolution, humign,data_filename,'humign')
1903	ENDIF
1904	!endspit
1905	ENDIF !! NOT disable_fire
1906
1907	IF (printlev_loc>=3) WRITE (numout,*) ' slowproc_init done '
1908
1909	END SUBROUTINE slowproc_init
1910
1911	!! ================================================================================================================================
1912	!! SUBROUTINE : slowproc_clear
1913	!!
1914	!>\BRIEF Clear all variables related to slowproc and stomate modules
1915	!!
1916	!_ ================================================================================================================================
1917
1918	SUBROUTINE slowproc_clear
1919
1920	! 1 clear all the variables defined as common for the routines in slowproc
1921
1922	IF (ALLOCATED (clayfraction)) DEALLOCATE (clayfraction)
1923	IF (ALLOCATED (laimap)) DEALLOCATE (laimap)
1924	IF (ALLOCATED (frac_nobio_lastyear)) DEALLOCATE (frac_nobio_lastyear)
1925	IF (ALLOCATED (vegetnew_firstday)) DEALLOCATE (vegetnew_firstday)
1926	IF (ALLOCATED (veget_max_new)) DEALLOCATE (veget_max_new)
1927	IF (ALLOCATED (frac_nobio_new)) DEALLOCATE (frac_nobio_new)
1928	IF ( ALLOCATED (soilclass_default)) DEALLOCATE (soilclass_default)
1929	!spitfire
1930	IF (ALLOCATED(m_lightn)) DEALLOCATE (m_lightn)
1931	IF (ALLOCATED(glccNetLCC)) DEALLOCATE (glccNetLCC)
1932	IF (ALLOCATED(glccSecondShift)) DEALLOCATE (glccSecondShift)
1933	IF (ALLOCATED(glccPrimaryShift)) DEALLOCATE (glccPrimaryShift)
1934	IF (ALLOCATED(harvest_matrix)) DEALLOCATE (harvest_matrix)
1935	IF (ALLOCATED(harvest_biomass)) DEALLOCATE (harvest_biomass)
1936	IF (ALLOCATED(bound_spa)) DEALLOCATE (bound_spa)
1937	IF (ALLOCATED(newvegfrac)) DEALLOCATE (newvegfrac)
1938	IF (ALLOCATED(proxy_anidens)) DEALLOCATE (proxy_anidens)
1939	IF (ALLOCATED(popd)) DEALLOCATE (popd)
1940	IF (ALLOCATED(humign)) DEALLOCATE (humign)
1941	IF (ALLOCATED(m_observed_ba)) DEALLOCATE (m_observed_ba)
1942	IF (ALLOCATED(m_cf_coarse)) DEALLOCATE (m_cf_coarse)
1943	IF (ALLOCATED(m_cf_fine)) DEALLOCATE (m_cf_fine)
1944	IF (ALLOCATED(m_ratio)) DEALLOCATE (m_ratio)
1945	IF (ALLOCATED(m_ratio_flag)) DEALLOCATE (m_ratio_flag)
1946	!endspit
1947	!gmjc fc_grazing
1948	IF (ALLOCATED (fc_grazing)) DEALLOCATE (fc_grazing)
1949	!end gmjc
1950	! 2. Clear all the variables in stomate
1951
1952	CALL stomate_clear
1953	!
1954	END SUBROUTINE slowproc_clear
1955
1956	!! ================================================================================================================================
1957	!! SUBROUTINE : slowproc_derivvar
1958	!!
1959	!>\BRIEF Initializes variables related to the
1960	!! parameters to be assimilated, the maximum water on vegetation, the vegetation height,
1961	!! and the fraction of soil covered by dead leaves and the vegetation height
1962	!!
1963	!! DESCRIPTION : (definitions, functional, design, flags):
1964	!! (1) Initialization of the variables relevant for the assimilation parameters
1965	!! (2) Intialization of the fraction of soil covered by dead leaves
1966	!! (3) Initialization of the Vegetation height per PFT
1967	!! (3) Initialization the maximum water on vegetation for interception with a particular treatement of the PFT no.1
1968	!!
1969	!! RECENT CHANGE(S): None
1970	!!
1971	!! MAIN OUTPUT VARIABLE(S): ::qsintmax, ::deadleaf_cover, ::assim_param, ::height
1972	!!
1973	!! REFERENCE(S) : None
1974	!!
1975	!! FLOWCHART : None
1976	!! \n
1977	!_ ================================================================================================================================
1978
1979	SUBROUTINE slowproc_derivvar (kjpindex, veget, lai, &
1980	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
1981
1982	!! INTERFACE DESCRIPTION
1983
1984	!! 0.1 Input scalar and fields
1985	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only
1986	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
1987	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! PFT leaf area index (m^{2} m^{-2})
1988
1989	!! 0.2. Output scalar and fields
1990	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water on vegetation for interception(mm)
1991	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover !! fraction of soil covered by dead leaves (unitless)
1992	REAL(r_std), DIMENSION (kjpindex,nvm,npco2), INTENT (out) :: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
1993	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height of the vegetation or surface in general ??? (m)
1994	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! growth temperature (°C)
1995	!
1996	!! 0.3 Local declaration
1997	INTEGER(i_std) :: jv !! Local indices
1998	!_ ================================================================================================================================
1999
2000	!
2001	! 1. Initialize (why here ??) the variables revelant for the assimilation parameters
2002	!
2003	DO jv = 1, nvm
2004	assim_param(:,jv,ivcmax) = vcmax_fix(jv)
2005	ENDDO
2006
2007	!
2008	! 2. Intialize the fraction of soil covered by dead leaves
2009	!
2010	deadleaf_cover(:) = zero
2011
2012	!
2013	! 3. Initialize the Vegetation height per PFT
2014	!
2015	DO jv = 1, nvm
2016	height(:,jv) = height_presc(jv)
2017	ENDDO
2018	!
2019	! 4. Initialize the maximum water on vegetation for interception
2020	!
2021	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
2022
2023	! Added by Nathalie - July 2006
2024	! Initialize the case of the PFT no.1 to zero
2025	qsintmax(:,1) = zero
2026
2027	temp_growth(:)=25.
2028
2029	END SUBROUTINE slowproc_derivvar
2030
2031
2032	!! ================================================================================================================================
2033	!! SUBROUTINE : slowproc_mean
2034	!!
2035	!>\BRIEF Accumulates field_in over a period of dt_tot.
2036	!! Has to be called at every time step (dt).
2037	!! Mean value is calculated if ldmean=.TRUE.
2038	!! field_mean must be initialized outside of this routine!
2039	!!
2040	!! DESCRIPTION : (definitions, functional, design, flags):
2041	!! (1) AcumAcuumlm
2042	!!
2043	!! RECENT CHANGE(S): None
2044	!!
2045	!! MAIN OUTPUT VARIABLE(S): ::field_main
2046	!!
2047	!! REFERENCE(S) : None
2048	!!
2049	!! FLOWCHART : None
2050	!! \n
2051	!_ ================================================================================================================================
2052
2053	SUBROUTINE slowproc_mean (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_mean)
2054
2055	!
2056	!! 0 declarations
2057
2058	!! 0.1 input scalar and variables
2059	INTEGER(i_std), INTENT(in) :: npts !! Domain size- terrestrial pixels only
2060	INTEGER(i_std), INTENT(in) :: n_dim2 !! Number of PFTs
2061	REAL(r_std), INTENT(in) :: dt_tot !! Time step of stomate (in days). The period over which the accumulation or the mean is computed
2062	REAL(r_std), INTENT(in) :: dt !! Time step in days
2063	LOGICAL, INTENT(in) :: ldmean !! Flag to calculate the mean after the accumulation ???
2064	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_in !! Daily field
2065
2066	!! 0.3 Modified field; The computed sum or mean field over dt_tot time period depending on the flag ldmean
2067	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_mean !! Accumulated field at dt_tot time period or mean field over dt_tot
2068
2069
2070	!_ ================================================================================================================================
2071
2072	!
2073	! 1. Accumulation the field over dt_tot period
2074	!
2075	field_mean(:,:) = field_mean(:,:) + field_in(:,:) * dt
2076
2077	!
2078	! 2. If the flag ldmean set, the mean field is computed over dt_tot period
2079	!
2080	IF (ldmean) THEN
2081	field_mean(:,:) = field_mean(:,:) / dt_tot
2082	ENDIF
2083
2084	END SUBROUTINE slowproc_mean
2085
2086
2087
2088	!! ================================================================================================================================
2089	!! SUBROUTINE : slowproc_long
2090	!!
2091	!>\BRIEF Calculates a temporally smoothed field (field_long) from
2092	!! instantaneous input fields.Time constant tau determines the strength of the smoothing.
2093	!! For tau -> infinity??, field_long becomes the true mean value of field_inst
2094	!! (but the spinup becomes infinietly long, too).
2095	!! field_long must be initialized outside of this routine!
2096	!!
2097	!! DESCRIPTION : (definitions, functional, design, flags):
2098	!! (1) Testing the time coherence betwen the time step dt and the time tau over which
2099	!! the rescaled of the mean is performed
2100	!! (2) Computing the rescaled mean over tau period
2101	!! MAIN OUTPUT VARIABLE(S): field_long
2102	!!
2103	!! RECENT CHANGE(S): None
2104	!!
2105	!! MAIN OUTPUT VARIABLE(S): ::field_long
2106	!!
2107	!! REFERENCE(S) : None
2108	!!
2109	!! FLOWCHART : None
2110	!! \n
2111	!_ ================================================================================================================================
2112
2113	SUBROUTINE slowproc_long (npts, n_dim2, dt, tau, field_inst, field_long)
2114
2115	!
2116	! 0 declarations
2117	!
2118
2119	! 0.1 input scalar and fields
2120
2121	INTEGER(i_std), INTENT(in) :: npts !! Domain size- terrestrial pixels only
2122	INTEGER(i_std), INTENT(in) :: n_dim2 !! Second dimension of the fields, which represents the number of PFTs
2123	REAL(r_std), INTENT(in) :: dt !! Time step in days
2124	REAL(r_std), INTENT(in) :: tau !! Integration time constant (has to have same unit as dt!)
2125	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_inst !! Instantaneous field
2126
2127
2128	! 0.2 modified field
2129
2130	! Long-term field
2131	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_long !! Mean value of the instantaneous field rescaled at tau time period
2132
2133	!_ ================================================================================================================================
2134
2135	!
2136	! 1 test coherence of the time
2137
2138	IF ( ( tau .LT. dt ) .OR. ( dt .LE. zero ) .OR. ( tau .LE. zero ) ) THEN
2139	WRITE(numout,*) 'slowproc_long: Problem with time steps'
2140	WRITE(numout,*) 'dt=',dt
2141	WRITE(numout,*) 'tau=',tau
2142	ENDIF
2143
2144	!
2145	! 2 integration of the field over tau
2146
2147	field_long(:,:) = ( field_inst(:,:)dt + field_long(:,:)(tau-dt) ) / tau
2148
2149	END SUBROUTINE slowproc_long
2150
2151
2152	!! ================================================================================================================================
2153	!! SUBROUTINE : slowproc_veget
2154	!!
2155	!>\BRIEF Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
2156	!!
2157	!! DESCRIPTION : Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
2158	!! (1) Set veget_max and frac_nobio for fraction smaller than min_vegfrac.
2159	!! (2) Reset some variables in stomate for small fractions
2160	!! (3) Calculate veget
2161	!! (5) Calculate totfrac_nobio
2162	!! (6) Calculate soiltile
2163	!!
2164	!! RECENT CHANGE(S): None
2165	!!
2166	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, totfrac_nobio, veget_max, veget, soiltile
2167	!!
2168	!! REFERENCE(S) : None
2169	!!
2170	!! FLOWCHART : None
2171	!! \n
2172	!_ ================================================================================================================================
2173
2174	SUBROUTINE slowproc_veget (kjpindex, f_rot_sech, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile)
2175	!
2176	! 0. Declarations
2177	!
2178	! 0.1 Input variables
2179	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
2180	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! PFT leaf area index (m^{2} m^{-2})
2181
2182	! 0.2 Modified variables
2183	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
2184	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
2185
2186	! 0.3 Output variables
2187	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
2188	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio
2189	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
2190	LOGICAL,DIMENSION(kjpindex), INTENT(in) :: f_rot_sech !! whether a grid point is under rotation
2191
2192	! 0.4 Local scalar and varaiables
2193	INTEGER(i_std) :: ji, jv, jst !! indices
2194	REAL(r_std) :: SUMveg
2195
2196	!_ ================================================================================================================================
2197	IF (printlev_loc > 8) WRITE(numout,*) 'Entering slowproc_veget'
2198
2199	! 0. Normalize fractions of frac_nobio and veget_max smaller than min_vegfrac
2200	! This is due to precision issues (float 64). It might lead to negative values.
2201	! At some point, a floating point exception.
2202	DO ji = 1, kjpindex
2203	IF ( SUM(frac_nobio(ji,:)) .LT. min_vegfrac ) THEN
2204	frac_nobio(ji,:) = zero
2205	ENDIF
2206
2207	IF (.NOT. ok_dgvm) THEN
2208	DO jv = 1, nvm
2209	IF ( veget_max(ji,jv) .LT. min_vegfrac ) THEN
2210	veget_max(ji,jv) = zero
2211	ENDIF
2212	ENDDO
2213	END IF
2214
2215	!! Normalize to keep the sum equal 1.
2216	SUMveg = SUM(frac_nobio(ji,:))+SUM(veget_max(ji,:))
2217	frac_nobio(ji,:) = frac_nobio(ji,:)/SUMveg
2218	veget_max(ji,:) = veget_max(ji,:)/SUMveg
2219	ENDDO
2220
2221
2222	!! 2. Reset some variables in stomate for small fractions.
2223	IF (ok_stomate .AND. .NOT. ok_dgvm) CALL stomate_veget_update(kjpindex,veget_max,f_rot_sech)
2224
2225
2226	!! 3. Calculate veget
2227	!! If lai of a vegetation type (jv > 1) is small, increase soil part
2228	!! stomate-like calculation
2229	DO ji = 1, kjpindex
2230	veget(ji,1)=veget_max(ji,1)
2231	DO jv = 2, nvm
2232	veget(ji,jv) = veget_max(ji,jv) * ( un - exp( - lai(ji,jv) * ext_coeff_vegetfrac(jv) ) )
2233	ENDDO
2234	ENDDO
2235
2236
2237	!! 4. Calculate totfrac_nobio
2238	totfrac_nobio(:) = zero
2239	DO jv = 1, nnobio
2240	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
2241	ENDDO
2242
2243
2244	!! 5. Calculate soiltiles
2245	!! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
2246	!! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
2247	!! in the other modules to perform separated energy balances
2248	! The sum of all soiltiles makes one, and corresponds to the bio fraction
2249	! of the grid cell (called vegtot in hydrol)
2250	soiltile(:,:) = zero
2251	DO jv = 1, nvm
2252	jst = pref_soil_veg(jv)
2253	DO ji = 1, kjpindex
2254	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
2255	ENDDO
2256	ENDDO
2257	DO ji = 1, kjpindex
2258	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
2259	soiltile(ji,:)=soiltile(ji,:)/(1.-totfrac_nobio(ji))
2260	ENDIF
2261	ENDDO
2262
2263	END SUBROUTINE slowproc_veget
2264
2265
2266	!! ================================================================================================================================
2267	!! SUBROUTINE : slowproc_lai
2268	!!
2269	!>\BRIEF Do the interpolation of lai for the PFTs in case the laimap is not read
2270	!!
2271	!! DESCRIPTION : (definitions, functional, design, flags):
2272	!! (1) Interplation by using the mean value of laimin and laimax for the PFTs
2273	!! (2) Interpolation between laimax and laimin values by using the temporal
2274	!! variations
2275	!! (3) If problem occurs during the interpolation, the routine stops
2276	!!
2277	!! RECENT CHANGE(S): None
2278	!!
2279	!! MAIN OUTPUT VARIABLE(S): ::lai
2280	!!
2281	!! REFERENCE(S) : None
2282	!!
2283	!! FLOWCHART : None
2284	!! \n
2285	!_ ================================================================================================================================
2286
2287	SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,laimap)
2288	!
2289	! 0. Declarations
2290	!
2291	!! 0.1 Input variables
2292	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
2293	INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI
2294	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag !! Soil temperature (K) ???
2295	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
2296	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x an y of the grid (m) - surface area of the gridbox
2297	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: laimap !! map of lai read
2298
2299	!! 0.2 Output
2300	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! PFT leaf area index (m^{2} m^{-2})LAI
2301
2302	!! 0.4 Local
2303	INTEGER(i_std) :: ji,jv !! Local indices
2304	!_ ================================================================================================================================
2305
2306	!
2307	IF ( .NOT. read_lai ) THEN
2308
2309	lai(: ,1) = zero
2310	! On boucle sur 2,nvm au lieu de 1,nvm
2311	DO jv = 2,nvm
2312	SELECT CASE (type_of_lai(jv))
2313
2314	CASE ("mean ")
2315	!
2316	! 1. do the interpolation between laimax and laimin
2317	!
2318	lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv))
2319	!
2320	CASE ("inter")
2321	!
2322	! 2. do the interpolation between laimax and laimin
2323	!
2324	DO ji = 1,kjpindex
2325	lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv))
2326	ENDDO
2327	!
2328	CASE default
2329	!
2330	! 3. Problem
2331	!
2332	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
2333	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
2334	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=false','')
2335	END SELECT
2336
2337	ENDDO
2338	!
2339	ELSE
2340	lai(: ,1) = zero
2341	! On boucle sur 2,nvm au lieu de 1,nvm
2342	DO jv = 2,nvm
2343
2344	SELECT CASE (type_of_lai(jv))
2345
2346	CASE ("mean ")
2347	!
2348	! 1. force MAXVAL of laimap on lai on this PFT
2349	!
2350	DO ji = 1,kjpindex
2351	lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
2352	ENDDO
2353	!
2354	CASE ("inter")
2355	!
2356	! 2. do the interpolation between laimax and laimin
2357	!
2358	!
2359	! If January
2360	!
2361	IF (month_end .EQ. 1 ) THEN
2362	IF (day_end .LE. 15) THEN
2363	lai(:,jv) = laimap(:,jv,12)(1-(day_end+15)/30.) + laimap(:,jv,1)((day_end+15)/30.)
2364	ELSE
2365	lai(:,jv) = laimap(:,jv,1)(1-(day_end-15)/30.) + laimap(:,jv,2)((day_end-15)/30.)
2366	ENDIF
2367	!
2368	! If December
2369	!
2370	ELSE IF (month_end .EQ. 12) THEN
2371	IF (day_end .LE. 15) THEN
2372	lai(:,jv) = laimap(:,jv,11)(1-(day_end+15)/30.) + laimap(:,jv,12)((day_end+15)/30.)
2373	ELSE
2374	lai(:,jv) = laimap(:,jv,12)(1-(day_end-15)/30.) + laimap(:,jv,1)((day_end-15)/30.)
2375	ENDIF
2376	!
2377	! ELSE
2378	!
2379	ELSE
2380	IF (day_end .LE. 15) THEN
2381	lai(:,jv) = laimap(:,jv,month_end-1)(1-(day_end+15)/30.) + laimap(:,jv,month_end)((day_end+15)/30.)
2382	ELSE
2383	lai(:,jv) = laimap(:,jv,month_end)(1-(day_end-15)/30.) + laimap(:,jv,month_end+1)((day_end-15)/30.)
2384	ENDIF
2385	ENDIF
2386	!
2387	CASE default
2388	!
2389	! 3. Problem
2390	!
2391	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
2392	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
2393	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=true','')
2394	END SELECT
2395
2396	ENDDO
2397	ENDIF
2398
2399	END SUBROUTINE slowproc_lai
2400
2401	!! ================================================================================================================================
2402	!! SUBROUTINE : slowproc_interlai
2403	!!
2404	!>\BRIEF Interpolate the LAI map to the grid of the model
2405	!!
2406	!! DESCRIPTION : (definitions, functional, design, flags):
2407	!!
2408	!! RECENT CHANGE(S): None
2409	!!
2410	!! MAIN OUTPUT VARIABLE(S): ::laimap
2411	!!
2412	!! REFERENCE(S) : None
2413	!!
2414	!! FLOWCHART : None
2415	!! \n
2416	!_ ================================================================================================================================
2417
2418	SUBROUTINE slowproc_interlai(nbpt, lalo, resolution, neighbours, contfrac, laimap)
2419
2420	USE interpweight
2421
2422	IMPLICIT NONE
2423
2424	!
2425	!
2426	!
2427	! 0.1 INPUT
2428	!
2429	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
2430	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
2431	!! (beware of the order = 1 : latitude, 2 : longitude)
2432	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
2433	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2434	!! (1=North and then clockwise)
2435	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
2436	!
2437	! 0.2 OUTPUT
2438	!
2439	REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) !! lai read variable and re-dimensioned
2440	!
2441	! 0.3 LOCAL
2442	!
2443	CHARACTER(LEN=80) :: filename !! name of the LAI map read
2444	INTEGER(i_std) :: ib, ip, jp, it, jv
2445	REAL(r_std) :: lmax, lmin, ldelta
2446	LOGICAL :: renormelize_lai ! flag to force LAI renormelization
2447	INTEGER :: ier
2448
2449	REAL(r_std), DIMENSION(nbpt) :: alaimap !! availability of the lai interpolation
2450	INTEGER, DIMENSION(4) :: invardims
2451	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE :: lairefrac !! lai fractions re-dimensioned
2452	REAL(r_std), DIMENSION(:), ALLOCATABLE :: vmin, vmax !! min/max values to use for the
2453	!! renormalization
2454	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2455	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2456	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
2457	!! (variabletypevals(1) = -un, not used)
2458	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2459	!! 'XYKindTime': Input values are kinds
2460	!! of something with a temporal
2461	!! evolution on the dx*dy matrix'
2462	LOGICAL :: nonegative !! whether negative values should be removed
2463	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2464	!! 'nomask': no-mask is applied
2465	!! 'mbelow': take values below maskvals(1)
2466	!! 'mabove': take values above maskvals(1)
2467	!! 'msumrange': take values within 2 ranges;
2468	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2469	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2470	!! (normalized by maskvals(3))
2471	!! 'var': mask values are taken from a
2472	!! variable inside the file (>0)
2473	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2474	!! `maskingtype')
2475	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2476	!_ ================================================================================================================================
2477
2478	!
2479	!Config Key = LAI_FILE
2480	!Config Desc = Name of file from which the vegetation map is to be read
2481	!Config If = LAI_MAP
2482	!Config Def = lai2D.nc
2483	!Config Help = The name of the file to be opened to read the LAI
2484	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
2485	!Config map which is derived from a Nicolas VIOVY one.
2486	!Config Units = [FILE]
2487	!
2488	filename = 'lai2D.nc'
2489	CALL getin_p('LAI_FILE',filename)
2490
2491	variablename = 'LAI'
2492
2493	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_interlai: Read and interpolate " &
2494	// TRIM(filename) //" for variable " //TRIM(variablename)
2495
2496	! invardims: shape of variable in input file to interpolate
2497	invardims = interpweight_get_var4dims_file(filename, variablename)
2498	! Check coherence of dimensions read from the file
2499	IF (invardims(4) /= 12) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of time dimension in input file for lai','','')
2500	IF (invardims(3) /= nvm) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of PFT dimension in input file for lai','','')
2501
2502	ALLOCATE(vmin(nvm),stat=ier)
2503	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmin','','')
2504
2505	ALLOCATE(vmax(nvm), STAT=ier)
2506	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmax','','')
2507
2508	ALLOCATE(lairefrac(nbpt,nvm,invardims(4)), STAT=ier)
2509	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable lairefrac','','')
2510
2511	! Assigning values to vmin, vmax
2512	vmin = un
2513	vmax = nvm*un
2514
2515	variabletypevals = -un
2516
2517	!! Variables for interpweight
2518	! Type of calculation of cell fractions
2519	fractype = 'default'
2520	! Name of the longitude and latitude in the input file
2521	lonname = 'longitude'
2522	latname = 'latitude'
2523	! Should negative values be set to zero from input file?
2524	nonegative = .TRUE.
2525	! Type of mask to apply to the input data (see header for more details)
2526	maskingtype = 'mbelow'
2527	! Values to use for the masking
2528	maskvals = (/ 20., undef_sechiba, undef_sechiba /)
2529	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2530	namemaskvar = ''
2531
2532	CALL interpweight_4D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
2533	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2534	maskvals, namemaskvar, nvm, invardims(4), -1, fractype, &
2535	-1., -1., lairefrac, alaimap)
2536
2537	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_interlai after interpweight_4D'
2538
2539	!
2540	!
2541	!Config Key = RENORM_LAI
2542	!Config Desc = flag to force LAI renormelization
2543	!Config If = LAI_MAP
2544	!Config Def = n
2545	!Config Help = If true, the laimap will be renormalize between llaimin and llaimax parameters.
2546	!Config Units = [FLAG]
2547	!
2548	renormelize_lai = .FALSE.
2549	CALL getin_p('RENORM_LAI',renormelize_lai)
2550
2551	!
2552	laimap(:,:,:) = zero
2553	!
2554	IF (printlev_loc >= 5) THEN
2555	WRITE(numout,*)' slowproc_interlai before starting loop nbpt:', nbpt
2556	END IF
2557
2558	! Assigning the right values and giving a value where information was not found
2559	DO ib=1,nbpt
2560	IF (alaimap(ib) < 0.) THEN
2561	DO jv=1,nvm
2562	laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux
2563	ENDDO
2564	ELSE
2565	DO jv=1, nvm
2566	DO it=1, invardims(4)
2567	laimap(ib,jv,it) = lairefrac(ib,jv,it)
2568	ENDDO
2569	ENDDO
2570	END IF
2571	ENDDO
2572	!
2573	! Normelize the read LAI by the values SECHIBA is used to
2574	!
2575	IF ( renormelize_lai ) THEN
2576	DO ib=1,nbpt
2577	DO jv=1, nvm
2578	lmax = MAXVAL(laimap(ib,jv,:))
2579	lmin = MINVAL(laimap(ib,jv,:))
2580	ldelta = lmax-lmin
2581	IF ( ldelta < min_sechiba) THEN
2582	! LAI constante ... keep it constant
2583	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)+(llaimax(jv)+llaimin(jv))/deux
2584	ELSE
2585	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)/(lmax-lmin)*(llaimax(jv)-llaimin(jv))+llaimin(jv)
2586	ENDIF
2587	ENDDO
2588	ENDDO
2589	ENDIF
2590
2591	! Write diagnostics
2592	CALL xios_orchidee_send_field("alaimap",alaimap)
2593
2594	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_interlai ended'
2595
2596	END SUBROUTINE slowproc_interlai
2597
2598	!! ================================================================================================================================
2599	!! SUBROUTINE : slowproc_readvegetmax
2600	!!
2601	!>\BRIEF Interpolate a vegetation map (by pft)
2602	!!
2603	!! DESCRIPTION : (definitions, functional, design, flags):
2604	!!
2605	!! RECENT CHANGE(S): None
2606	!!
2607	!! MAIN OUTPUT VARIABLE(S):
2608	!!
2609	!! REFERENCE(S) : None
2610	!!
2611	!! FLOWCHART : None
2612	!! \n
2613	!_ ================================================================================================================================
2614
2615	SUBROUTINE slowproc_readvegetmax(nbpt, lalo, neighbours, resolution, contfrac, &
2616	veget_last, veget_next, frac_nobio_next, veget_year, init)
2617
2618	USE interpweight
2619
2620	IMPLICIT NONE
2621
2622	!
2623	!
2624	!
2625	! 0.1 INPUT
2626	!
2627	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs
2628	!! to be interpolated
2629	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
2630	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
2631	!! (1=North and then clockwise)
2632	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
2633	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
2634	!
2635	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last !! old max vegetfrac
2636	INTEGER(i_std), INTENT(in) :: veget_year !! first year for landuse (0 == NO TIME AXIS)
2637	LOGICAL, INTENT(in) :: init !! initialisation
2638	!! In case of dgvm == FALSE, whatever its value,
2639	!! all PFT fractions will be updated.
2640	!
2641	! 0.2 OUTPUT
2642	!
2643	REAL(r_std), DIMENSION(nbpt,nvmap), INTENT(out) :: veget_next !! new max vegetfrac
2644	REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next !! new fraction of the mesh which is
2645	!! covered by ice, lakes, ...
2646
2647	!
2648	! 0.3 LOCAL
2649	!
2650	!
2651	CHARACTER(LEN=80) :: filename
2652	INTEGER(i_std) :: ib, inobio, jv
2653	REAL(r_std) :: sumf, err, norm
2654	!
2655	! for DGVM case :
2656	REAL(r_std) :: sum_veg ! sum of vegets
2657	REAL(r_std) :: sum_nobio ! sum of nobios
2658	REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets
2659	REAL(r_std) :: rapport ! (S-B) / (S-A)
2660	LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones)
2661	! e.g. in case of DGVM and not init (optional parameter)
2662	REAL(r_std), DIMENSION(nbpt,nvm) :: vegetrefrac !! veget fractions re-dimensioned
2663	REAL(r_std), DIMENSION(nbpt) :: aveget !! Availability of the soilcol interpolation
2664	REAL(r_std), DIMENSION(nvm) :: vmin, vmax !! min/max values to use for the renormalization
2665	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2666	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2667	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
2668	!! (variabletypevals(1) = -un, not used)
2669	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2670	!! 'XYKindTime': Input values are kinds
2671	!! of something with a temporal
2672	!! evolution on the dx*dy matrix'
2673	LOGICAL :: nonegative !! whether negative values should be removed
2674	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2675	!! 'nomask': no-mask is applied
2676	!! 'mbelow': take values below maskvals(1)
2677	!! 'mabove': take values above maskvals(1)
2678	!! 'msumrange': take values within 2 ranges;
2679	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2680	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2681	!! (normalized by maskvals(3))
2682	!! 'var': mask values are taken from a
2683	!! variable inside the file (>0)
2684	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2685	!! `maskingtype')
2686	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2687	CHARACTER(LEN=250) :: msg
2688
2689	!_ ================================================================================================================================
2690
2691	IF (printlev_loc >= 5) PRINT *,' In slowproc_readvegetmax'
2692
2693	!
2694	!Config Key = VEGETATION_FILE
2695	!Config Desc = Name of file from which the vegetation map is to be read
2696	!Config If = MAP_PFT_FORMAT
2697	!Config Def = PFTmap.nc
2698	!Config Help = The name of the file to be opened to read a vegetation
2699	!Config map (in pft) is to be given here.
2700	!Config Units = [FILE]
2701	!
2702	filename = 'PFTmap.nc'
2703	CALL getin_p('VEGETATION_FILE',filename)
2704	variablename = 'maxvegetfrac'
2705
2706	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readvegetmax: Read and interpolate " &
2707	// TRIM(filename) // " for variable " // TRIM(variablename)
2708
2709	! Assigning values to vmin, vmax
2710	vmin = 1
2711	! chaoyue: be careful of this when do merging with trunk. This has to be
2712	! persistently diverged from trunk because trunk does not have, neither uses
2713	! the concept of age calsses. Here nvmap is the number of metaclasses.
2714	vmax = nvmap*1._r_std
2715
2716	variabletypevals = -un
2717
2718	!! Variables for interpweight
2719	! Type of calculation of cell fractions
2720	fractype = 'default'
2721	! Name of the longitude and latitude in the input file
2722	lonname = 'lon'
2723	latname = 'lat'
2724	! Should negative values be set to zero from input file?
2725	nonegative = .FALSE.
2726	! Type of mask to apply to the input data (see header for more details)
2727	maskingtype = 'msumrange'
2728	! Values to use for the masking
2729	maskvals = (/ 1.-1.e-7, 0., 2. /)
2730	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2731	namemaskvar = ''
2732
2733	! persistent divergence with trunk. nvmap rather than nvm should be used.
2734	! see the comment above.
2735	CALL interpweight_3D(nbpt, nvmap, variabletypevals, lalo, resolution, neighbours, &
2736	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2737	maskvals, namemaskvar, nvmap, 0, veget_year, fractype, &
2738	-1., -1., vegetrefrac, aveget)
2739	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after interpeeight_3D'
2740
2741	!
2742	! Compute the logical for partial (only anthropic) PTFs update
2743	IF (ok_dgvm .AND. .NOT. init) THEN
2744	partial_update= .TRUE.
2745	ELSE
2746	partial_update=.FALSE.
2747	END IF
2748
2749	IF (printlev_loc >= 5) THEN
2750	WRITE(numout,*)' slowproc_readvegetmax before updating loop nbpt:', nbpt
2751	END IF
2752
2753	IF ( .NOT. partial_update ) THEN
2754	veget_next(:,:)=zero
2755
2756	IF (printlev_loc >=3 .AND. ANY(aveget < min_sechiba)) THEN
2757	WRITE(numout,*) 'Some grid cells on the model grid did not have any points on the source grid.'
2758	IF (init) THEN
2759	WRITE(numout,*) 'Initialization with full fraction of bare soil are done for the below grid cells.'
2760	ELSE
2761	WRITE(numout,*) 'Old values are kept for the below grid cells.'
2762	ENDIF
2763	WRITE(numout,*) 'List of grid cells (ib, lat, lon):'
2764	END IF
2765
2766	DO ib = 1, nbpt
2767	! vegetrefrac is already normalized to sum equal one for each grid cell
2768	veget_next(ib,:) = vegetrefrac(ib,:)
2769
2770	IF (aveget(ib) < min_sechiba) THEN
2771	IF (printlev_loc >=3) WRITE(numout,*) ib,lalo(ib,1),lalo(ib,2)
2772	IF (init) THEN
2773	veget_next(ib,1) = un
2774	veget_next(ib,2:nvmap) = zero
2775	ELSE
2776	veget_next(ib,:) = veget_last(ib,:)
2777	ENDIF
2778	ENDIF
2779	ENDDO
2780	! `partial_update` is TRUE, the natural PFT fraction will be copied from
2781	! `veget_last` to `veget_next`, only the anthropogenic PFT are updated
2782	! from the input veget_max map. This is used only when DGVM is activated.
2783	!print *,"second time within slowproc_readvegetmax"
2784	!print *,"veget_last",veget_last
2785	!print *,"veget_next",veget_next
2786	!print *,"partial_update",partial_update
2787	!print *,"vegmap",vegmap
2788	ELSE
2789	! Partial update
2790	DO ib = 1, nbpt
2791	IF (aveget(ib) > min_sechiba) THEN
2792	! For the case with properly interpolated grid cells (aveget>0)
2793
2794	! last veget for this point
2795	sum_veg=SUM(veget_last(ib,:))
2796	!
2797	! If the DGVM is activated, only anthropic PFTs are utpdated, the others are copied from previous time-step
2798	veget_next(ib,:) = veget_last(ib,:)
2799
2800	DO jv = 2, nvmap
2801	IF ( .NOT. natural(jv) .OR. pasture(jv)) THEN
2802	veget_next(ib,jv) = vegetrefrac(ib,jv)
2803	ENDIF
2804	ENDDO
2805
2806	sumvAnthro_old = zero
2807	sumvAnthro = zero
2808	DO jv = 2, nvmap
2809	IF ( .NOT. natural(jv) .OR. pasture(jv)) THEN
2810	sumvAnthro = sumvAnthro + veget_next(ib,jv)
2811	sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv)
2812	ENDIF
2813	ENDDO
2814
2815	IF ( sumvAnthro_old < sumvAnthro ) THEN
2816	! Increase of non natural vegetations (increase of agriculture)
2817	! The proportion of natural PFT's must be preserved
2818	! ie the sum of vegets is preserved
2819	! and natural PFT / (sum of veget - sum of antropic veget)
2820	! is preserved.
2821	rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old )
2822	DO jv = 1, nvmap
2823	IF ( natural(jv) .AND. .NOT. pasture(jv)) THEN
2824	veget_next(ib,jv) = veget_last(ib,jv) * rapport
2825	ENDIF
2826	ENDDO
2827	ELSE
2828	! Increase of natural vegetations (decrease of agriculture)
2829	! The decrease of agriculture is replaced by bare soil. The DGVM will
2830	! re-introduce natural PFT's.
2831	DO jv = 1, nvmap
2832	IF ( natural(jv) .AND. .NOT. pasture(jv)) THEN
2833	veget_next(ib,jv) = veget_last(ib,jv)
2834	ENDIF
2835	ENDDO
2836	veget_next(ib,1) = veget_next(ib,1) + sumvAnthro_old - sumvAnthro
2837	ENDIF
2838
2839	! test
2840	IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > 10*EPSILON(un) ) THEN
2841	WRITE(numout,*) 'slowproc_readvegetmax _______'
2842	msg = " No conservation of sum of veget for point "
2843	WRITE(numout,*) TRIM(msg), ib, ",(", lalo(ib,1),",", lalo(ib,2), ")"
2844	WRITE(numout,*) " last sum of veget ", sum_veg, " new sum of veget ", &
2845	SUM(veget_next(ib,:)), " error : ", SUM(veget_next(ib,:))-sum_veg
2846	WRITE(numout,*) " Anthropic modifications : last ",sumvAnthro_old," new ",sumvAnthro
2847	CALL ipslerr_p (3,'slowproc_readvegetmax', &
2848	& 'No conservation of sum of veget_next', &
2849	& "The sum of veget_next is different after reading Land Use map.", &
2850	& '(verify the dgvm case model.)')
2851	ENDIF
2852	ELSE
2853	! For the case when there was a propblem with the interpolation, aveget < min_sechiba
2854	WRITE(numout,*) 'slowproc_readvegetmax _______'
2855	WRITE(numout,*) " No land point in the map for point ", ib, ",(", lalo(ib,1), ",", &
2856	lalo(ib,2),")"
2857	CALL ipslerr_p (2,'slowproc_readvegetmax', &
2858	& 'Problem with vegetation file for Land Use.', &
2859	& "No land point in the map for point", &
2860	& '(verify your land use file.)')
2861	veget_next(ib,:) = veget_last(ib,:)
2862	ENDIF
2863
2864	ENDDO
2865	ENDIF
2866
2867	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after updating'
2868	!
2869	frac_nobio_next (:,:) = un
2870	!
2871	! Work only for one nnobio !! (ie ice)
2872	DO inobio=1,nnobio
2873	DO jv=1,nvmap
2874	DO ib = 1, nbpt
2875	frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv)
2876	ENDDO
2877	ENDDO
2878	ENDDO
2879
2880	DO ib = 1, nbpt
2881	sum_veg = SUM(veget_next(ib,:))
2882	sum_nobio = SUM(frac_nobio_next(ib,:))
2883	IF (sum_nobio < 0.) THEN
2884	frac_nobio_next(ib,:) = zero
2885	veget_next(ib,1) = veget_next(ib,1) + sum_nobio
2886	sum_veg = SUM(veget_next(ib,:))
2887	ENDIF
2888	sumf = sum_veg + sum_nobio
2889	IF (sumf > min_sechiba) THEN
2890	veget_next(ib,:) = veget_next(ib,:) / sumf
2891	frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf
2892	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2893	err=norm-un
2894	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ",ib, &
2895	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un, sumf",err,sumf
2896	IF (abs(err) > -EPSILON(un)) THEN
2897	IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN
2898	frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err
2899	ELSE
2900	veget_next(ib,1) = veget_next(ib,1) - err
2901	ENDIF
2902	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2903	err=norm-un
2904	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ", ib, &
2905	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un",err
2906	IF (abs(err) > EPSILON(un)) THEN
2907	WRITE(numout,*) ' slowproc_readvegetmax _______'
2908	WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2909	WRITE(numout,*) " err(sum-1.) = ",abs(err)
2910	CALL ipslerr_p (2,'slowproc_readvegetmax', &
2911	& 'Problem with sum vegetation + sum fracnobio for Land Use.', &
2912	& "sum not equal to 1.", &
2913	& '(verify your land use file.)')
2914	aveget(ib) = -0.6
2915	ENDIF
2916	ENDIF
2917	ELSE
2918	! sumf < min_sechiba
2919	WRITE(numout,*) ' slowproc_readvegetmax _______'
2920	WRITE(numout,*)" No vegetation nor frac_nobio for point ", ib, ",(", lalo(ib,1), ",", &
2921	lalo(ib,2),")"
2922	WRITE(numout,*)" Replaced by bare_soil !! "
2923	veget_next(ib,1) = un
2924	veget_next(ib,2:nvmap) = zero
2925	frac_nobio_next(ib,:) = zero
2926	ENDIF
2927	ENDDO
2928
2929	! Write diagnostics
2930	CALL xios_orchidee_send_field("aveget",aveget)
2931
2932	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_readvegetmax ended'
2933
2934	END SUBROUTINE slowproc_readvegetmax
2935
2936	!! ================================================================================================================================
2937	!! SUBROUTINE : slowproc_interpol
2938	!!
2939	!>\BRIEF Interpolate the IGBP vegetation map to the grid of the model
2940	!!
2941	!! DESCRIPTION : (definitions, functional, design, flags):
2942	!!
2943	!! RECENT CHANGE(S): None
2944	!!
2945	!! MAIN OUTPUT VARIABLE(S):
2946	!!
2947	!! REFERENCE(S) : None
2948	!!
2949	!! FLOWCHART : None
2950	!! \n
2951	!_ ================================================================================================================================
2952
2953	SUBROUTINE slowproc_interpol(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio)
2954
2955	USE interpweight
2956
2957	IMPLICIT NONE
2958
2959	!
2960	!
2961	!
2962	! 0.1 INPUT
2963	!
2964	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
2965	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order!)
2966	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2967	!! (1=North and then clockwise)
2968	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
2969	REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
2970	!
2971	! 0.2 OUTPUT
2972	!
2973	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) !! Vegetation fractions
2974	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) !! Fraction of the mesh which is covered by ice, lakes, ...
2975	!
2976	! 0.3 LOCAL
2977	!
2978	INTEGER(i_std), PARAMETER :: nolson = 94 !! Number of Olson classes
2979	REAL(r_std) :: resollon, resollat !! resolution of the longitudes and the latitudes
2980	!! in the input data which it is in a Goode compressed projection
2981	!
2982	!
2983	CHARACTER(LEN=80) :: filename
2984	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, vid
2985	REAL(r_std), DIMENSION(1) :: lev
2986	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap
2987	REAL(r_std) :: vegcorr(nolson,nvm)
2988	REAL(r_std) :: nobiocorr(nolson,nnobio)
2989	REAL(r_std) :: sumf
2990	INTEGER(i_std) :: jv, inear
2991	INTEGER :: ALLOC_ERR
2992	INTEGER :: Ndimslonlat !! Number of dimensions of lon/lat
2993	CHARACTER(LEN=1) :: dimlLS
2994	INTEGER :: dim1Dlonlat !! Length of 1D longitudes, latitudes
2995	INTEGER, DIMENSION(2) :: invardims2D
2996	REAL(r_std), DIMENSION(nbpt,nolson) :: vegetrefrac !! vegegt fractions re-dimensioned
2997	REAL(r_std), DIMENSION(nbpt) :: aveget5k !! Availability of the interpolation
2998	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: aveget5k_glob !! Availability of the interpolation
2999	REAL(r_std) :: vmin, vmax !! min/max values to use for the
3000	!! renormalization
3001	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3002	CHARACTER(LEN=80) :: lonname, latname !! lonm lat names in input file
3003	REAL(r_std), DIMENSION(nolson) :: variabletypevals !! Values for all the types of the variable
3004	!! (variabletypevals(1) = -un, not used)
3005	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3006	!! 'XYKindTime': Input values are kinds
3007	!! of something with a temporal
3008	!! evolution on the dx*dy matrix'
3009	LOGICAL :: nonegative !! whether negative values should be removed
3010	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3011	!! 'nomask': no-mask is applied
3012	!! 'mbelow': take values below maskvals(1)
3013	!! 'mabove': take values above maskvals(1)
3014	!! 'msumrange': take values within 2 ranges;
3015	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3016	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3017	!! (normalized by maskvals(3))
3018	!! 'var': mask values are taken from a
3019	!! variable inside the file (>0)
3020	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3021	!! `maskingtype')
3022	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3023	LOGICAL :: foundnegvals !! whether negative aveget5k values
3024	!! where found
3025	CHARACTER(LEN=250) :: msg
3026	INTEGER :: rcode
3027
3028	!_ ================================================================================================================================
3029	variablename = 'vegetation_map'
3030
3031	CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3032	!Config Key = VEGETATION_FILE
3033	!Config Desc = Name of file from which the vegetation map is to be read
3034	!Config If = NOT(IMPOSE_VEG) and NOT(MAP_PFT_FORMAT)
3035	!Config Def = carteveg5km.nc
3036	!Config Help = The name of the file to be opened to read the vegetation
3037	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
3038	!Config map which is derived from the IGBP one. We assume that we have
3039	!Config a classification in 87 types. This is Olson modified by Viovy.
3040	!Config Units = [FILE]
3041	!
3042	filename = 'carteveg5km.nc'
3043	CALL getin_p('VEGETATION_FILE',filename)
3044
3045	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_interpol: Read and interpolate " &
3046	// TRIM(filename) // " for variable " // TRIM(variablename)
3047
3048	! Assigning values to vmin, vmax
3049	vmin = un
3050	vmax = nolson*un
3051	!
3052	!
3053	ALLOC_ERR=-1
3054	variabletypevals = -un
3055
3056	!! Variables for interpweight
3057	! Type of calculation of cell fractions
3058	fractype = 'default'
3059	! Name of the longitude and latitude in the input file
3060	lonname = 'longitude'
3061	latname = 'latitude'
3062	! Should negative values be set to zero from input file?
3063	nonegative = .FALSE.
3064	! Type of mask to apply to the input data (see header for more details)
3065	maskingtype = 'mabove'
3066	! Values to use for the masking
3067	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3068	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
3069	namemaskvar = ''
3070	! Meridional and zonal resolutions of the input data [m]
3071	resollon = 5000.*un
3072	resollat = 5000.*un
3073
3074	CALL interpweight_1D(nbpt, nolson, variabletypevals, lalo, resolution, neighbours, &
3075	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3076	maskvals, namemaskvar, 0, 0, -1, fractype, &
3077	resollon, resollat, vegetrefrac, aveget5k)
3078	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_interpol after interpweight_1D'
3079
3080
3081	!
3082	! Some assumptions on the vegetation file. This information should be
3083	! be computed or read from the file.
3084	! It is the resolution in meters of the grid of the vegetation file.
3085	!
3086	!
3087	! Now we know how many points of which Olson type from the fine grid fall
3088	! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3089	!
3090	! vegetrefrac is already normalized in subroutine interpweight_1D
3091	!
3092	! now finally calculate coarse vegetation map
3093	! Find which model vegetation corresponds to each Olson type
3094	!
3095	veget(:,:) = zero
3096	frac_nobio(:,:) = zero
3097
3098	DO vid = 1, nolson
3099	DO jv = 1, nvm
3100	veget(:,jv) = veget(:,jv) + vegetrefrac(:,vid) * vegcorr(vid,jv)
3101	ENDDO
3102
3103	DO jv = 1, nnobio
3104	frac_nobio(:,jv) = frac_nobio(:,jv) + vegetrefrac(:,vid) * nobiocorr(vid,jv)
3105	ENDDO
3106	ENDDO
3107
3108	IF (printlev_loc >= 5) THEN
3109	WRITE(numout,*)' slowproc_interpol before starting loop nbpt:', nbpt
3110	END IF
3111
3112	! Getting input longitude and latitude matrices for looking nearest cell
3113	! Looking on the global grid if there are points without interpolated values
3114	IF (is_root_prc) THEN
3115	ALLOCATE(aveget5k_glob(iim_g*jjm_g))
3116	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_interpol','Problem in allocation of variable aveget5k_glo','','')
3117	ELSE
3118	ALLOCATE (aveget5k_glob(1))
3119	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_interpol','Problem in allocation of variable aveget5k_glo','','')
3120	ENDIF
3121	CALL gather(aveget5k,aveget5k_glob)
3122
3123	IF (is_root_prc) THEN
3124	foundnegvals = ANY(aveget5k_glob .lt. zero)
3125	END IF
3126	CAll bcast(foundnegvals)
3127
3128	IF (foundnegvals) THEN
3129	! lon, lat matrices of the input data have to be recupered...
3130	WRITE(numout,*) ' Looking for nearest point on the 5 km map'
3131	IF (is_root_prc) THEN
3132	CALL flininfo(filename, dim1Dlonlat, jml, lml, tml, fid)
3133	!Ndimslonlat = interpweight_get_varNdims_file(filename, TRIM(lonname))
3134	END IF
3135	!CALL bcast(Ndimslonlat)
3136	CALL bcast(dim1Dlonlat)
3137	Ndimslonlat = 1
3138	IF (Ndimslonlat ==1) THEN
3139	ALLOCATE(lon_ful(dim1Dlonlat), STAT=ALLOC_ERR)
3140	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_interpol','Problem in allocation of variable lon_ful','','')
3141	ALLOCATE(lat_ful(dim1Dlonlat), STAT=ALLOC_ERR)
3142	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_interpol','Problem in allocation of variable lat_ful','','')
3143	ALLOCATE(vegmap(dim1Dlonlat), STAT=ALLOC_ERR)
3144	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_interpol','Problem in allocation of variable vegmap','','')
3145	IF (is_root_prc) THEN
3146	CALL flinget(fid, TRIM(lonname), dim1Dlonlat, 0, 0, 0, 1, 1, lon_ful)
3147	CALL flinget(fid, TRIM(latname), dim1Dlonlat, 0, 0, 0, 1, 1, lat_ful)
3148	CALL flinget(fid, TRIM(variablename), dim1Dlonlat, 0, 0, 0, 1, 1, vegmap)
3149	CALL flinclo(fid)
3150	END IF
3151	ELSE
3152	WRITE(dimlLS,'(A1)')Ndimslonlat
3153	msg = "Problem in rank of '" // TRIM(lonname) // "': " // dimlLS // " Not ready !!"
3154	CALL ipslerr_p(3,'slowproc_interpol',TRIM(msg),'','')
3155	END IF
3156	END IF
3157
3158	CALL bcast(lon_ful)
3159	CALL bcast(lat_ful)
3160	CALL bcast(vegmap)
3161
3162	DEALLOCATE(aveget5k_glob)
3163
3164	!
3165	! Clean up the point of the map
3166	!
3167	DO ib = 1, nbpt
3168	!
3169	! Let us see if all points found something in the 5km map !
3170	!
3171	IF ( aveget5k(ib) .EQ. -1 ) THEN
3172	!
3173	! Now we need to handle some exceptions
3174	!
3175	IF ( lalo(ib,1) .LT. -56.0) THEN
3176	! Antartica
3177	frac_nobio(ib,:) = zero
3178	frac_nobio(ib,iice) = un
3179	veget(ib,:) = zero
3180	! aveget5k(ib) = -1.2
3181	ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3182	! Artica
3183	frac_nobio(ib,:) = zero
3184	frac_nobio(ib,iice) = un
3185	veget(ib,:) = zero
3186	! aveget5k(ib) = -1.2
3187	ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3188	! Greenland
3189	frac_nobio(ib,:) = zero
3190	frac_nobio(ib,iice) = un
3191	veget(ib,:) = zero
3192	! aveget5k(ib) = -1.2
3193	ELSE
3194	WRITE(numout,*) ' slowproc_interpol _______'
3195	WRITE(numout,*) ' PROBLEM, no point in the 5km map found for this grid box',ib
3196	WRITE(numout,*) ' Longitude range : ', lalo(ib,2)
3197	WRITE(numout,*) ' Latitude range : ', lalo(ib,1)
3198
3199	CALL slowproc_nearest (dim1Dlonlat, lon_ful, lat_ful, &
3200	lalo(ib,2), lalo(ib,1), inear)
3201	WRITE(numout,*) ' Coordinates of the nearest point:', &
3202	lon_ful(inear),lat_ful(inear)
3203
3204	DO jv = 1, nvm
3205	veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3206	ENDDO
3207
3208	DO jv = 1, nnobio
3209	frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3210	ENDDO
3211	ENDIF
3212	ENDIF
3213	!
3214	!
3215	! Limit the smallest vegetation fraction to 0.5%
3216	!
3217	DO vid = 1, nvm
3218	IF ( veget(ib,vid) .LT. min_vegfrac ) THEN
3219	veget(ib,vid) = zero
3220	ENDIF
3221	ENDDO
3222
3223	sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3224	frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3225	veget(ib,:) = veget(ib,:)/sumf
3226	ENDDO
3227
3228	IF (ALLOCATED(vegmap)) DEALLOCATE(vegmap)
3229	IF (ALLOCATED(lon_ful)) DEALLOCATE(lon_ful)
3230	IF (ALLOCATED(lat_ful)) DEALLOCATE(lat_ful)
3231
3232	! Write diagnostics
3233	CALL xios_orchidee_send_field("aveget5k",aveget5k)
3234
3235	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_interpol ended'
3236
3237	END SUBROUTINE slowproc_interpol
3238
3239	!! ================================================================================================================================
3240	!! SUBROUTINE : slowproc_interpol_g
3241	!!
3242	!>\BRIEF Interpolate the IGBP vegetation map to the grid of the model
3243	!!
3244	!! DESCRIPTION : (definitions, functional, design, flags):
3245	!!
3246	!! RECENT CHANGE(S): None
3247	!!
3248	!! MAIN OUTPUT VARIABLE(S): ::veget, ::frac_nobio
3249	!!
3250	!! REFERENCE(S) : None
3251	!!
3252	!! FLOWCHART : None
3253	!! \n
3254	!_ ================================================================================================================================
3255
3256	SUBROUTINE slowproc_interpol_g(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
3257	!
3258	!
3259	!
3260	! 0.1 INPUT
3261	!
3262	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
3263	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
3264	!! (beware of the order : 1=latitude ; 2=longitude)
3265	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point
3266	!! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
3267	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
3268	REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
3269	!
3270	! 0.2 OUTPUT
3271	!
3272	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) !! Vegetation fractions
3273	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) !! Fraction of the mesh which is covered by ice, lakes, ...
3274	!
3275	LOGICAL :: ok_interpol !! optionnal return of aggregate_vec
3276	!
3277	! 0.3 LOCAL
3278	!
3279	INTEGER(i_std), PARAMETER :: nolson = 94 !! Number of Olson classes
3280	!
3281	!
3282	CHARACTER(LEN=80) :: filename !!vegetation map filename
3283	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, vid
3284	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful, vegmap !! for 5km vegetation map
3285	!! latitude vector, longitude vector, and
3286	!! value of Olson's classes for each location
3287	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the model grid ???
3288	!! cf src_global/interpol_help.f90, line 377, called "areaoverlap"
3289	INTEGER(i_std),ALLOCATABLE, DIMENSION(:,:) :: sub_index !! the indexes from the grid boxes from the data that go
3290	!! into the model's boxes
3291	!! cf src_global/interpol_help.f90,line 300, called "ip"
3292	REAL(r_std), DIMENSION(nbpt,nolson) :: n_origveg !! number of points of each Olson type from the fine grid
3293	!! in each box of the (coarse) model grid
3294	REAL(r_std), DIMENSION(nbpt) :: n_found !! total number of different Olson types found in each
3295	!! box of the (coarse) model grid
3296	REAL(r_std), DIMENSION(nbpt,nolson) :: frac_origveg !! fraction of each Olson type in each box of the (coarse) model grid
3297	REAL(r_std) :: vegcorr(nolson,nvm) !! correspondance table between Olson and the following SECHIBA Classes.
3298	!! vegcorr(i,:)+nobiocorr(i,:) = 1. for all i
3299	!! see each class in src_parameters/constantes_veg.f90
3300
3301	REAL(r_std) :: nobiocorr(nolson,nnobio) !! non-biospheric surface typesi
3302	CHARACTER(LEN=40) :: callsign !! Allows to specify which variable is beeing treated
3303	REAL(r_std) :: sumf, resol_lon, resol_lat !! sumf = sum veget + sum nobio
3304	!! resol_lon, resol_lat reolution in meters of the grid of the vegetation file
3305	INTEGER(i_std) :: idi, jv, inear, nbvmax !! idi : counter for nbvmax, see below
3306	!! jv : counter for nvm, number of PFT
3307	!! inear : location of the point of vegmap, which is the closest from the modelled point
3308	!! nbvmax : number of maximum vegetation map points in the GCM grid
3309	INTEGER(i_std) :: nix, njx
3310	!
3311	INTEGER :: ALLOC_ERR !! location of the eventual missing value in vegmap
3312
3313	!_ ================================================================================================================================
3314	!
3315	n_origveg(:,:) = zero
3316	n_found(:) = zero
3317	!
3318	CALL get_vegcorr (nolson,vegcorr,nobiocorr)
3319	!
3320	!Config Key = VEGETATION_FILE
3321	!Config Desc = Name of file from which the vegetation map is to be read
3322	!Config If = NOT(IMPOSE_VEG) and NOT(MAP_PFT_FORMAT)
3323	!Config Def = carteveg5km.nc
3324	!Config Help = The name of the file to be opened to read the vegetation
3325	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
3326	!Config map which is derived from the IGBP one. We assume that we have
3327	!Config a classification in 87 types. This is Olson modified by Viovy.
3328	!Config Units = [FILE]
3329	!
3330	filename = 'carteveg5km.nc'
3331	CALL getin_p('VEGETATION_FILE',filename) ! GETIN_P !!
3332	!
3333	CALL flininfo(filename, iml, jml, lml, tml, fid)
3334	!
3335	! see IOIPSL/src/flincom.f90, line 665
3336	! fid : File ID
3337	!- iml \| These 4 variables give the size of the variables
3338	!- jml \| to be read. It will be verified that the variables
3339	!- lml \| fits in there.
3340	!- tml
3341	! iml, jml : horizontal size of the grid, lml = vertical size
3342	! tml : size of time axis
3343
3344	! TL : pourquoi 2 variables pour la taille horizontale ? cf
3345	! IOIPSL/src/flincom.f90 , line 160
3346
3347	ALLOC_ERR=-1
3348	ALLOCATE(lat_ful(iml), STAT=ALLOC_ERR)
3349	IF (ALLOC_ERR/=0) THEN
3350	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
3351	STOP
3352	ENDIF
3353	ALLOC_ERR=-1
3354	ALLOCATE(lon_ful(iml), STAT=ALLOC_ERR)
3355	IF (ALLOC_ERR/=0) THEN
3356	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
3357	STOP
3358	ENDIF
3359	ALLOC_ERR=-1
3360	ALLOCATE(vegmap(iml), STAT=ALLOC_ERR)
3361	IF (ALLOC_ERR/=0) THEN
3362	WRITE(numout,*) "ERROR IN ALLOCATION of vegmap : ",ALLOC_ERR
3363	STOP
3364	ENDIF
3365	!
3366	WRITE(numout,*) 'Reading the OLSON type vegetation file'
3367	!
3368	CALL flinget(fid, 'longitude', iml, jml, lml, tml, 1, 1, lon_ful)
3369	CALL flinget(fid, 'latitude', iml, jml, lml, tml, 1, 1, lat_ful)
3370	CALL flinget(fid, 'vegetation_map', iml, jml, lml, tml, 1, 1, vegmap)
3371	!
3372	WRITE(numout,*) 'File name : ', filename
3373	WRITE(numout,*) 'Min and max vegetation numbers : ', MINVAL(vegmap), MAXVAL(vegmap)
3374	!
3375	CALL flinclo(fid)
3376	!
3377	IF (MAXVAL(vegmap) .LT. nolson) THEN
3378	WRITE(numout,*) 'WARNING -- WARNING'
3379	WRITE(numout,*) 'The vegetation map has too few vegetation types.'
3380	WRITE(numout,*) 'If you are lucky it will work but please check'
3381	ELSE IF ( MAXVAL(vegmap) .GT. nolson) THEN
3382	WRITE(numout,*) 'More vegetation types in file than the code can'
3383	WRITE(numout,*) 'deal with.: ', MAXVAL(vegmap), nolson
3384	STOP 'slowproc_interpol'
3385	ENDIF
3386	!
3387	! Some assumptions on the vegetation file. This information should be
3388	! be computed or read from the file.
3389	! It is the reolution in meters of the grid of the vegetation file.
3390	!
3391
3392	!TL : CODE EN DUR ?????
3393	resol_lon = 5000.
3394	resol_lat = 5000.
3395	!
3396	!
3397	! The number of maximum vegetation map points in the GCM grid is estimated.
3398	! Some margin is taken.
3399	!
3400	nix=INT(MAXVAL(resolution_g(:,1)*2)/resol_lon)+1
3401	njx=INT(MAXVAL(resolution_g(:,2)*2)/resol_lon)+1
3402	nbvmax = nix*njx
3403	!
3404	! No need to broadcast as this routine is only called on root_proc
3405	!
3406	callsign="Vegetation map"
3407	!
3408	ok_interpol = .FALSE.
3409	DO WHILE ( .NOT. ok_interpol )
3410	WRITE(numout,*) "Projection arrays for ",callsign," : "
3411	WRITE(numout,*) "nbvmax = ",nbvmax
3412	!
3413	ALLOC_ERR=-1
3414	ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
3415	IF (ALLOC_ERR/=0) THEN
3416	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
3417	STOP
3418	ENDIF
3419	sub_index(:,:)=0
3420	ALLOC_ERR=-1
3421	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
3422	IF (ALLOC_ERR/=0) THEN
3423	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
3424	STOP
3425	ENDIF
3426	sub_area(:,:)=zero
3427	!
3428	CALL aggregate_p (nbpt, lalo, neighbours, resolution, contfrac, &
3429	& iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, &
3430	& nbvmax, sub_index, sub_area, ok_interpol)
3431	!
3432	! Defined as aggregate_2d or aggregate_vec in src_global/interpol_help.f90, depending
3433	! on the dimensions (2D region or vector)i.
3434	! This routing will get for each point of the coarse grid the
3435	! indexes of the finer grid and the area of overlap.
3436	! This routine is designed for a fine grid which is regular in lat/lon.
3437
3438	IF ( .NOT. ok_interpol ) THEN
3439	DEALLOCATE(sub_area)
3440	DEALLOCATE(sub_index)
3441	!
3442	nbvmax = nbvmax * 2
3443	ELSE
3444	!
3445	DO ib = 1, nbpt
3446	DO idi=1, nbvmax
3447	! Leave the do loop if all sub areas are treated, sub_area <= 0
3448	IF ( sub_area(ib,idi) <= zero ) EXIT
3449
3450	ip = sub_index(ib,idi)
3451	n_origveg(ib,NINT(vegmap(ip))) = n_origveg(ib,NINT(vegmap(ip))) + sub_area(ib,idi)
3452	n_found(ib) = n_found(ib) + sub_area(ib,idi)
3453	ENDDO
3454	ENDDO
3455	!
3456	ENDIF
3457	ENDDO
3458	!
3459	! Now we know how many points of which Olson type from the fine grid fall
3460	! into each box of the (coarse) model grid: n_origveg(nbpt,nolson)
3461	!
3462	!
3463	! determine fraction of Olson vegetation type in each box of the coarse grid
3464	!
3465	DO vid = 1, nolson
3466	WHERE ( n_found(:) .GT. 0 )
3467	frac_origveg(:,vid) = n_origveg(:,vid) / n_found(:)
3468	ELSEWHERE
3469	frac_origveg(:,vid) = zero
3470	ENDWHERE
3471	ENDDO
3472	!
3473	! now finally calculate coarse vegetation map
3474	! Find which model vegetation corresponds to each Olson type
3475	!
3476	veget(:,:) = zero
3477	frac_nobio(:,:) = zero
3478	!
3479	DO vid = 1, nolson
3480	!
3481	DO jv = 1, nvm
3482	veget(:,jv) = veget(:,jv) + frac_origveg(:,vid) * vegcorr(vid,jv)
3483	ENDDO
3484	!
3485	DO jv = 1, nnobio
3486	frac_nobio(:,jv) = frac_nobio(:,jv) + frac_origveg(:,vid) * nobiocorr(vid,jv)
3487	ENDDO
3488	!
3489	ENDDO
3490	!
3491	WRITE(numout,*) 'slowproc_interpol : Interpolation Done'
3492	!
3493	! Clean up the point of the map
3494	!
3495	DO ib = 1, nbpt
3496	!
3497	! Let us see if all points found something in the 5km map !
3498	!
3499	IF ( n_found(ib) .EQ. 0 ) THEN
3500	!
3501	! Now we need to handle some exceptions
3502	!
3503	IF ( lalo(ib,1) .LT. -56.0) THEN
3504	! Antartica
3505	frac_nobio(ib,:) = zero
3506	frac_nobio(ib,iice) = un
3507	veget(ib,:) = zero
3508	!
3509	ELSE IF ( lalo(ib,1) .GT. 70.0) THEN
3510	! Artica
3511	frac_nobio(ib,:) = zero
3512	frac_nobio(ib,iice) = un
3513	veget(ib,:) = zero
3514	!
3515	ELSE IF ( lalo(ib,1) .GT. 55.0 .AND. lalo(ib,2) .GT. -65.0 .AND. lalo(ib,2) .LT. -20.0) THEN
3516	! Greenland
3517	frac_nobio(ib,:) = zero
3518	frac_nobio(ib,iice) = un
3519	veget(ib,:) = zero
3520	!
3521	ELSE
3522	!
3523	WRITE(numout,*) 'PROBLEM, no point in the 5km map found for this grid box',ib
3524	WRITE(numout,*) 'Longitude range : ', lalo(ib,2)
3525	WRITE(numout,*) 'Latitude range : ', lalo(ib,1)
3526	!
3527	WRITE(numout,*) 'Looking for nearest point on the 5 km map'
3528	CALL slowproc_nearest (iml, lon_ful, lat_ful, &
3529	lalo(ib,2), lalo(ib,1), inear)
3530	WRITE(numout,*) 'Coordinates of the nearest point:', &
3531	lon_ful(inear),lat_ful(inear)
3532	!
3533	DO jv = 1, nvm
3534	veget(ib,jv) = vegcorr(NINT(vegmap(inear)),jv)
3535	ENDDO
3536	!
3537	DO jv = 1, nnobio
3538	frac_nobio(ib,jv) = nobiocorr(NINT(vegmap(inear)),jv)
3539	ENDDO
3540	!
3541	ENDIF
3542	!
3543	ENDIF
3544	!
3545	!
3546	! Limit the smallest vegetation fraction to 0.5%
3547	!
3548	DO vid = 1, nvm
3549	IF ( veget(ib,vid) .LT. min_vegfrac ) THEN ! min_vegfrac=0.001 in constantes_veg.f90
3550	veget(ib,vid) = zero
3551	ENDIF
3552	ENDDO
3553	!
3554	sumf = SUM(frac_nobio(ib,:))+SUM(veget(ib,:))
3555	frac_nobio(ib,:) = frac_nobio(ib,:)/sumf
3556	veget(ib,:) = veget(ib,:)/sumf
3557	!
3558	!
3559	ENDDO
3560	!
3561	DEALLOCATE(vegmap)
3562	DEALLOCATE(lat_ful, lon_ful)
3563	DEALLOCATE(sub_index)
3564	DEALLOCATE(sub_area)
3565
3566	!
3567	RETURN
3568	!
3569	END SUBROUTINE slowproc_interpol_g
3570
3571
3572	!! ================================================================================================================================
3573	!! SUBROUTINE : slowproc_nearest
3574	!!
3575	!>\BRIEF looks for nearest grid point on the fine map
3576	!!
3577	!! DESCRIPTION : (definitions, functional, design, flags):
3578	!!
3579	!! RECENT CHANGE(S): None
3580	!!
3581	!! MAIN OUTPUT VARIABLE(S): ::inear
3582	!!
3583	!! REFERENCE(S) : None
3584	!!
3585	!! FLOWCHART : None
3586	!! \n
3587	!_ ================================================================================================================================
3588
3589	SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
3590
3591	!! INTERFACE DESCRIPTION
3592
3593	!! 0.1 input variables
3594
3595	INTEGER(i_std), INTENT(in) :: iml !! size of the vector
3596	REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5 !! longitude and latitude vector, for the 5km vegmap
3597	REAL(r_std), INTENT(in) :: lonmod, latmod !! longitude and latitude modelled
3598
3599	!! 0.2 output variables
3600
3601	INTEGER(i_std), INTENT(out) :: inear !! location of the grid point from the 5km vegmap grid
3602	!! closest from the modelled grid point
3603
3604	!! 0.4 Local variables
3605
3606	REAL(r_std) :: pa, p
3607	REAL(r_std) :: coscolat, sincolat
3608	REAL(r_std) :: cospa, sinpa
3609	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
3610	INTEGER(i_std) :: i
3611	INTEGER(i_std), DIMENSION(1) :: ineartab
3612	INTEGER :: ALLOC_ERR
3613
3614	!_ ================================================================================================================================
3615
3616	ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
3617	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_nearest','Error in allocation for cosang','','')
3618
3619	pa = pi/2.0 - latmod*pi/180.0 ! dist. between north pole and the point a
3620	!! COLATITUDE, in radian
3621	cospa = COS(pa)
3622	sinpa = SIN(pa)
3623
3624	DO i = 1, iml
3625
3626	sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 ) !! sinus of the colatitude
3627	coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 ) !! cosinus of the colatitude
3628
3629	p = (lonmod-lon5(i))*pi/180.0 !! angle between a & b (between their meridian)in radians
3630
3631	!! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
3632	cosang(i) = cospacoscolat + sinpasincolat*COS(p) !! TL : cosang is maximum when angle is at minimal value
3633	!! orthodromic distance between 2 points : cosang = cosinus (arc(AB)/R), with
3634	!R = Earth radius, then max(cosang) = max(cos(arc(AB)/R)), reached when arc(AB)/R is minimal, when
3635	! arc(AB) is minimal, thus when point B (corresponding grid point from LAI MAP) is the nearest from
3636	! modelled A point
3637	ENDDO
3638
3639	ineartab = MAXLOC( cosang(:) )
3640	inear = ineartab(1)
3641
3642	DEALLOCATE(cosang)
3643	END SUBROUTINE slowproc_nearest
3644
3645	!! ================================================================================================================================
3646	!! SUBROUTINE : slowproc_soilt
3647	!!
3648	!>\BRIEF Interpolate the Zobler soil type map
3649	!!
3650	!! DESCRIPTION : (definitions, functional, design, flags):
3651	!!
3652	!! RECENT CHANGE(S): None
3653	!!
3654	!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction
3655	!!
3656	!! REFERENCE(S) : None
3657	!!
3658	!! FLOWCHART : None
3659	!! \n
3660	!_ ================================================================================================================================
3661	SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
3662
3663	USE interpweight
3664
3665	IMPLICIT NONE
3666	!
3667	!
3668	! This subroutine should read the Zobler map and interpolate to the model grid. The method
3669	! is to get fraction of the three main soiltypes for each grid box.
3670	! The soil fraction are going to be put into the array soiltype in the following order :
3671	! coarse, medium and fine.
3672	!
3673	!
3674	!! 0.1 INPUT
3675	!
3676	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
3677	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
3678	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
3679	!! (1=North and then clockwise)
3680	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
3681	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
3682	!
3683	! 0.2 OUTPUT
3684	!
3685	REAL(r_std), INTENT(out) :: soilclass(nbpt, nscm) !! Soil type map to be created from the Zobler map
3686	REAL(r_std), INTENT(out) :: clayfraction(nbpt) !! The fraction of clay as used by STOMATE
3687	!
3688	!
3689	! 0.3 LOCAL
3690	!
3691	CHARACTER(LEN=80) :: filename
3692	INTEGER(i_std) :: ib, ilf, nbexp, i
3693	INTEGER(i_std) :: fopt !! Nb of pts from the texture map within one ORCHIDEE grid-cell
3694	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: soiltext, soiltext2
3695	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
3696	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
3697	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt !! Texture the different points from the input texture map
3698	!! in one ORCHIDEE grid cell (unitless)
3699	!
3700	! Number of texture classes in Zobler
3701	!
3702	INTEGER(i_std), PARAMETER :: nzobler = 7
3703	REAL(r_std),ALLOCATABLE :: textfrac_table(:,:)
3704	!
3705	INTEGER :: ALLOC_ERR
3706	INTEGER :: ntextinfile !! number of soil textures in the in the file
3707	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: textrefrac !! text fractions re-dimensioned
3708	REAL(r_std), DIMENSION(nbpt) :: atext !! Availability of the text interpolation
3709	REAL(r_std) :: vmin, vmax !! min/max values to use for the
3710
3711	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3712	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in input file
3713	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
3714	!! (variabletypevals(1) = -un, not used)
3715	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3716	!! 'XYKindTime': Input values are kinds
3717	!! of something with a temporal
3718	!! evolution on the dx*dy matrix'
3719	LOGICAL :: nonegative !! whether negative values should be removed
3720	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3721	!! 'nomask': no-mask is applied
3722	!! 'mbelow': take values below maskvals(1)
3723	!! 'mabove': take values above maskvals(1)
3724	!! 'msumrange': take values within 2 ranges;
3725	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3726	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3727	!! (normalized by maskvals(3))
3728	!! 'var': mask values are taken from a
3729	!! variable inside the file (>0)
3730	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3731	!! `maskingtype')
3732	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3733	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
3734	REAL(r_std) :: sgn !! sum of fractions excluding glaciers and ocean
3735	!_ ================================================================================================================================
3736
3737	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt'
3738	!
3739	! Needs to be a configurable variable
3740	!
3741	!
3742	!Config Key = SOILCLASS_FILE
3743	!Config Desc = Name of file from which soil types are read
3744	!Config Def = soils_param.nc
3745	!Config If = NOT(IMPOSE_VEG)
3746	!Config Help = The name of the file to be opened to read the soil types.
3747	!Config The data from this file is then interpolated to the grid of
3748	!Config of the model. The aim is to get fractions for sand loam and
3749	!Config clay in each grid box. This information is used for soil hydrology
3750	!Config and respiration.
3751	!Config Units = [FILE]
3752	!
3753	filename = 'soils_param.nc'
3754	CALL getin_p('SOILCLASS_FILE',filename)
3755
3756	variablename = 'soiltext'
3757
3758	!! Variables for interpweight
3759	! Type of calculation of cell fractions
3760	fractype = 'default'
3761	! Name of the longitude and latitude in the input file
3762	lonname = 'nav_lon'
3763	latname = 'nav_lat'
3764
3765	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
3766	// TRIM(filename) // " for variable " // TRIM(variablename)
3767
3768	IF ( TRIM(soil_classif) /= 'none' ) THEN
3769
3770	! Define a variable for the number of soil textures in the input file
3771	SELECTCASE(soil_classif)
3772	CASE('zobler')
3773	ntextinfile=nzobler
3774	CASE('usda')
3775	ntextinfile=nscm
3776	CASE DEFAULT
3777	WRITE(numout,*) 'slowproc_soilt:'
3778	WRITE(numout,*) ' A non supported soil type classification has been chosen'
3779	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
3780	ENDSELECT
3781
3782	ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR)
3783	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',&
3784	'','')
3785
3786	! Assigning values to vmin, vmax
3787	vmin = un
3788	vmax = ntextinfile*un
3789
3790	ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR)
3791	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','')
3792	variabletypevals = -un
3793
3794	!! Variables for interpweight
3795	! Should negative values be set to zero from input file?
3796	nonegative = .FALSE.
3797	! Type of mask to apply to the input data (see header for more details)
3798	maskingtype = 'mabove'
3799	! Values to use for the masking
3800	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3801	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
3802	namemaskvar = ''
3803
3804	CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours, &
3805	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3806	maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext)
3807
3808	ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR)
3809	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','')
3810	ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR)
3811	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','')
3812
3813	IF (printlev_loc >= 5) THEN
3814	WRITE(numout,*)' slowproc_soilt after interpweight_2D'
3815	WRITE(numout,*)' slowproc_soilt before starting loop nbpt:', nbpt
3816	WRITE(numout,*)" slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..."
3817	END IF
3818	ELSE
3819	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt using default values all points are propertly ' // &
3820	'interpolated atext = 1. everywhere!'
3821	atext = 1.
3822	END IF
3823
3824	nbexp = 0
3825	SELECTCASE(soil_classif)
3826	CASE('none')
3827	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
3828	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3829	DO ib=1, nbpt
3830	soilclass(ib,:) = soilclass_default_fao
3831	clayfraction(ib) = clayfraction_default
3832	ENDDO
3833	CASE('zobler')
3834	!
3835	soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
3836	!
3837	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
3838	!
3839	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
3840	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3841	CALL get_soilcorr_zobler (nzobler, textfrac_table)
3842	!
3843	!
3844	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt after getting table of textures'
3845	DO ib =1, nbpt
3846	soilclass(ib,:) = zero
3847	clayfraction(ib) = zero
3848	!
3849	! vecpos: List of positions where textures were not zero
3850	! vecpos(1): number of not null textures found
3851	vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq')
3852	fopt = vecpos(1)
3853
3854	IF ( fopt .EQ. 0 ) THEN
3855	! No points were found for current grid box, use default values
3856	nbexp = nbexp + 1
3857	soilclass(ib,:) = soilclass_default(:)
3858	clayfraction(ib) = clayfraction_default
3859	ELSE
3860	IF (fopt == nzobler) THEN
3861	! All textures are not zero
3862	solt=(/(i,i=1,nzobler)/)
3863	ELSE
3864	DO ilf = 1,fopt
3865	solt(ilf) = vecpos(ilf+1)
3866	END DO
3867	END IF
3868	!
3869	! Compute the fraction of each textural class
3870	!
3871	sgn = 0.
3872	DO ilf = 1,fopt
3873	!
3874	! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE
3875	! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine)
3876	! For type 6 = glacier, default values are set and it is also taken into account during the normalization
3877	! of the fractions (done in interpweight_2D)
3878	! Note that type 0 corresponds to ocean but it is already removed using the mask above.
3879	!
3880	IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. &
3881	(solt(ilf) .NE. 6) ) THEN
3882	SELECT CASE(solt(ilf))
3883	CASE(1)
3884	soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf))
3885	CASE(2)
3886	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3887	CASE(3)
3888	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3889	CASE(4)
3890	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3891	CASE(5)
3892	soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf))
3893	CASE(7)
3894	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3895	CASE DEFAULT
3896	WRITE(numout,*) 'We should not be here, an impossible case appeared'
3897	CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','')
3898	END SELECT
3899	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
3900	! over the zobler pixels composing the ORCHIDEE grid-cell
3901	clayfraction(ib) = clayfraction(ib) + &
3902	& textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf))
3903	! Sum the fractions which are not glaciers nor ocean
3904	sgn = sgn + textrefrac(ib,solt(ilf))
3905	ELSE
3906	IF (solt(ilf) .GT. nzobler) THEN
3907	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
3908	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','')
3909	ENDIF
3910	ENDIF
3911	ENDDO
3912
3913	IF ( sgn .LT. min_sechiba) THEN
3914	! Set default values if grid cells were only covered by glaciers or ocean
3915	! or if now information on the source grid was found.
3916	nbexp = nbexp + 1
3917	soilclass(ib,:) = soilclass_default(:)
3918	clayfraction(ib) = clayfraction_default
3919	ELSE
3920	! Normalize using the fraction of surface not including glaciers and ocean
3921	soilclass(ib,:) = soilclass(ib,:)/sgn
3922	clayfraction(ib) = clayfraction(ib)/sgn
3923	ENDIF
3924	ENDIF
3925	ENDDO
3926
3927	!
3928	CASE("fao")
3929	!
3930	soilclass_default=soilclass_default_fao
3931	!
3932	WRITE(numout,*) "Using a soilclass map with fao classification"
3933	!
3934	ALLOCATE(textfrac_table(nscm,ntext))
3935	!
3936	CALL get_soilcorr_zobler (nscm, textfrac_table)
3937	!
3938	DO ib =1, nbpt
3939	!
3940	! GO through the point we have found
3941	!
3942	!
3943	fopt = COUNT(sub_area(ib,:) > zero)
3944	!
3945	! Check that we found some points
3946	!
3947	soilclass(ib,:) = 0.0
3948	clayfraction(ib) = 0.0
3949	!
3950	IF ( fopt .EQ. 0) THEN
3951	nbexp = nbexp + 1
3952	soilclass(ib,:) = soilclass_default(:)
3953	clayfraction(ib) = clayfraction_default
3954	ELSE
3955	!
3956	DO ilf = 1,fopt
3957	solt(ilf) = soiltext(sub_index(ib,ilf,1),sub_index(ib,ilf,2))
3958	ENDDO
3959	!
3960	!
3961	! Compute the average bare soil albedo parameters
3962	!
3963	sgn = zero
3964	!
3965	DO ilf = 1,fopt
3966	!
3967	!
3968	!
3969	IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
3970	soilclass(ib,solt(ilf)) = soilclass(ib,solt(ilf)) + sub_area(ib,ilf)
3971	clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * sub_area(ib,ilf)
3972	sgn = sgn + sub_area(ib,ilf)
3973	ELSE
3974	IF (solt(ilf) .GT. nscm) THEN
3975	WRITE(,) 'The file contains a soil color class which is incompatible with this program'
3976	STOP 'slowproc_soilt'
3977	ENDIF
3978	ENDIF
3979	!
3980	ENDDO
3981	!
3982	! Normalize the surface
3983	!
3984	IF ( sgn .LT. min_sechiba) THEN
3985	nbexp = nbexp + 1
3986	soilclass(ib,:) = soilclass_default(:)
3987	clayfraction(ib) = clayfraction_default
3988	ELSE
3989	soilclass(ib,:) = soilclass(ib,:)/sgn
3990	clayfraction(ib) = clayfraction(ib)/sgn
3991	ENDIF
3992	!
3993	ENDIF
3994	!
3995	ENDDO
3996
3997	! The "USDA" case reads a map of the 12 USDA texture classes,
3998	! such as to assign the corresponding soil properties
3999	CASE("usda")
4000	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification"
4001
4002	soilclass_default=soilclass_default_usda
4003
4004	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
4005	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
4006
4007	CALL get_soilcorr_usda (nscm, textfrac_table)
4008
4009	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
4010	!
4011	DO ib =1, nbpt
4012	!
4013	! GO through the point we have found
4014	!
4015	!
4016	! Provide which textures were found
4017	! vecpos: List of positions where textures were not zero
4018	! vecpos(1): number of not null textures found
4019	vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq')
4020	fopt = vecpos(1)
4021
4022	!
4023	! Check that we found some points
4024	!
4025	soilclass(ib,:) = 0.0
4026	clayfraction(ib) = 0.0
4027
4028	IF ( fopt .EQ. 0) THEN
4029	! No points were found for current grid box, use default values
4030	IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib
4031	nbexp = nbexp + 1
4032	soilclass(ib,:) = soilclass_default
4033	clayfraction(ib) = clayfraction_default
4034	ELSE
4035	IF (fopt == nscm) THEN
4036	! All textures are not zero
4037	solt(:) = (/(i,i=1,nscm)/)
4038	ELSE
4039	DO ilf = 1,fopt
4040	solt(ilf) = vecpos(ilf+1)
4041	END DO
4042	END IF
4043	!
4044	!
4045	! Compute the fraction of each textural class
4046	!
4047	!
4048	DO ilf = 1,fopt
4049	IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
4050	soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf))
4051	clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * &
4052	textrefrac(ib,solt(ilf))
4053	ELSE
4054	IF (solt(ilf) .GT. nscm) THEN
4055	WRITE(,) 'The file contains a soil color class which is incompatible with this program'
4056	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','')
4057	ENDIF
4058	ENDIF
4059	!
4060	ENDDO
4061
4062	! Set default values if the surface in source file is too small
4063	IF ( atext(ib) .LT. min_sechiba) THEN
4064	nbexp = nbexp + 1
4065	soilclass(ib,:) = soilclass_default(:)
4066	clayfraction(ib) = clayfraction_default
4067	ENDIF
4068	ENDIF
4069
4070	ENDDO
4071
4072	IF (printlev_loc>=4) WRITE (numout,*) ' slowproc_soilt: End case usda'
4073
4074	CASE DEFAULT
4075	WRITE(numout,*) 'slowproc_soilt _______'
4076	WRITE(numout,*) ' A non supported soil type classification has been chosen'
4077	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
4078	ENDSELECT
4079	IF (printlev_loc >= 5 ) WRITE(numout,*)' slowproc_soilt end of type classification'
4080
4081	IF ( nbexp .GT. 0 ) THEN
4082	WRITE(numout,*) 'slowproc_soilt:'
4083	WRITE(numout,*) ' The interpolation of the bare soil albedo had ', nbexp
4084	WRITE(numout,*) ' points without data. This are either coastal points or ice covered land.'
4085	WRITE(numout,*) ' The problem was solved by using the default soil types.'
4086	ENDIF
4087
4088	IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals)
4089	IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac)
4090	IF (ALLOCATED(solt)) DEALLOCATE (solt)
4091	IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table)
4092
4093	! Write diagnostics
4094	CALL xios_orchidee_send_field("atext",atext)
4095
4096	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_soilt ended'
4097
4098	END SUBROUTINE slowproc_soilt
4099
4100	!! ================================================================================================================================
4101	!! SUBROUTINE : slowproc_slope
4102	!!
4103	!>\BRIEF Calculate mean slope coef in each model grid box from the slope map
4104	!!
4105	!! DESCRIPTION : (definitions, functional, design, flags):
4106	!!
4107	!! RECENT CHANGE(S): None
4108	!!
4109	!! MAIN OUTPUT VARIABLE(S): ::reinf_slope
4110	!!
4111	!! REFERENCE(S) : None
4112	!!
4113	!! FLOWCHART : None
4114	!! \n
4115	!_ ================================================================================================================================
4116
4117	SUBROUTINE slowproc_slope(nbpt, lalo, neighbours, resolution, contfrac, reinf_slope)
4118
4119	USE interpweight
4120
4121	IMPLICIT NONE
4122
4123	!
4124	!
4125	!
4126	! 0.1 INPUT
4127	!
4128	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
4129	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
4130	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point
4131	! (1=North and then clockwise)
4132	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
4133	REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Fraction of continent in the grid
4134	!
4135	! 0.2 OUTPUT
4136	!
4137	REAL(r_std), INTENT(out) :: reinf_slope(nbpt) ! slope coef
4138	!
4139	! 0.3 LOCAL
4140	!
4141	!
4142	REAL(r_std) :: slope_noreinf ! Slope above which runoff is maximum
4143	CHARACTER(LEN=80) :: filename
4144	REAL(r_std) :: vmin, vmax !! min/max values to use for the
4145	!! renormalization
4146	REAL(r_std), DIMENSION(nbpt) :: aslope !! slope availability
4147
4148	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
4149	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
4150	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
4151	!! 'XYKindTime': Input values are kinds
4152	!! of something with a temporal
4153	!! evolution on the dx*dy matrix'
4154	LOGICAL :: nonegative !! whether negative values should be removed
4155	CHARACTER(LEN=50) :: maskingtype !! Type of masking
4156	!! 'nomask': no-mask is applied
4157	!! 'mbelow': take values below maskvals(1)
4158	!! 'mabove': take values above maskvals(1)
4159	!! 'msumrange': take values within 2 ranges;
4160	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
4161	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
4162	!! (normalized by maskvals(3))
4163	!! 'var': mask values are taken from a
4164	!! variable inside the file (>0)
4165	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
4166	!! `maskingtype')
4167	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
4168
4169	!_ ================================================================================================================================
4170
4171	!
4172	!Config Key = SLOPE_NOREINF
4173	!Config Desc = See slope_noreinf above
4174	!Config If =
4175	!Config Def = 0.5
4176	!Config Help = The slope above which there is no reinfiltration
4177	!Config Units = [-]
4178	!
4179	slope_noreinf = 0.5
4180	!
4181	CALL getin_p('SLOPE_NOREINF',slope_noreinf)
4182	!
4183	!Config Key = TOPOGRAPHY_FILE
4184	!Config Desc = Name of file from which the topography map is to be read
4185	!Config If =
4186	!Config Def = cartepente2d_15min.nc
4187	!Config Help = The name of the file to be opened to read the orography
4188	!Config map is to be given here. Usualy SECHIBA runs with a 2'
4189	!Config map which is derived from the NGDC one.
4190	!Config Units = [FILE]
4191	!
4192	filename = 'cartepente2d_15min.nc'
4193	CALL getin_p('TOPOGRAPHY_FILE',filename)
4194
4195	variablename = 'pente'
4196	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_slope: Read and interpolate " &
4197	// TRIM(filename) // " for variable " // TRIM(variablename)
4198
4199	! For this case there are not types/categories. We have 'only' a continuos field
4200	! Assigning values to vmin, vmax
4201	vmin = 0.
4202	vmax = 9999.
4203
4204	!! Variables for interpweight
4205	! Type of calculation of cell fractions
4206	fractype = 'slopecalc'
4207	! Name of the longitude and latitude in the input file
4208	lonname = 'longitude'
4209	latname = 'latitude'
4210	! Should negative values be set to zero from input file?
4211	nonegative = .FALSE.
4212	! Type of mask to apply to the input data (see header for more details)
4213	maskingtype = 'mabove'
4214	! Values to use for the masking
4215	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
4216	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
4217	namemaskvar = ''
4218
4219	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
4220	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
4221	maskvals, namemaskvar, -1, fractype, slope_default, slope_noreinf, &
4222	reinf_slope, aslope)
4223	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_slope after interpweight_2Dcont'
4224
4225	! Write diagnostics
4226	CALL xios_orchidee_send_field("aslope",aslope)
4227
4228	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_slope ended'
4229
4230	END SUBROUTINE slowproc_slope
4231
4232	!! ================================================================================================================================
4233	!! SUBROUTINE : get_vegcorr
4234	!!
4235	!>\BRIEF The "get_vegcorr" routine defines the table of correspondence
4236	!! between the 94 Olson vegetation types and the 13 Plant Functional Types known
4237	!! by SECHIBA and STOMATE. Used by slowproc for the old interpolation.
4238	!!
4239	!!\DESCRIPTION : get_vegcorr is needed if you use the old_map carteveg5km.nc. \n
4240	!! Usually SECHIBA can run with a 5kmx5km map which is derived from the IGBP one. \n
4241	!! We assume that we have a classification in 94 types. This is Olson one modified by Nicolas Viovy.\n
4242	!! ORCHIDEE has to convert the Olson vegetation types into PFTs for the run (interpolation step).\n
4243	!! Each Olson matches to a combination of fractions of one or several PFTs.\n
4244	!! This routine uses the same process for the non-biospheric map (not used).\n
4245	!!
4246	!! RECENT CHANGE(S): None
4247	!!
4248	!! MAIN OUTPUT VARIABLE(S): ::vegcorr, ::nobiocorr.
4249	!!
4250	!! REFERENCE(S) :
4251	!! - Olson, J.S., J.A. Watts, and L.J. Allison., 1983.
4252	!! "Carbon in Live Vegetation of Major World Ecosystems."
4253	!! Report ORNL-5862. Oak Ridge National Laboratory, Oak Ridge, Tennessee.
4254	!! - Olson, J.S., J.A. Watts, and L.J. Allison., 1985.
4255	!! "Major World Ecosystem Complexes Ranked by Carbon in Live Vegetation: A Database."
4256	!! NDP-017. Carbon Dioxide Information Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee.
4257	!!
4258	!! FLOWCHART : None
4259	!! \n
4260	!_ ================================================================================================================================
4261
4262	SUBROUTINE get_vegcorr (nolson,vegcorr,nobiocorr)
4263
4264	IMPLICIT NONE
4265
4266	!! 0. Variables and parameters declaration
4267
4268	INTEGER(i_std),PARAMETER :: nolson94 = 94 !! Number of Olson vegetation types (unitless)
4269	INTEGER(i_std),PARAMETER :: nvm13 = 13 !! Number of PFTS of ORCHIDEE (unitless)
4270
4271	!! 0.1 Input variables
4272
4273	INTEGER(i_std),INTENT(in) :: nolson !! Number of Olson vegetation types (unitless)
4274
4275	!! 0.2 Output variables
4276
4277	REAL(r_std),DIMENSION(nolson,nvm),INTENT(out) :: vegcorr !! Correspondence array between Olson types and PFTS
4278	!! (0-1, unitless)
4279	REAL(r_std),DIMENSION(nolson,nnobio),INTENT(out) :: nobiocorr !! Correspondence array between non-vegetation types and nobio
4280	!! types (lake,etc..) (0-1, unitless)
4281
4282	!! 0.4 Local variable
4283
4284	INTEGER(i_std) :: ib !! Indice (unitless)
4285
4286	!_ ================================================================================================================================
4287
4288	!-
4289	! 0. Check consistency
4290	!-
4291	IF (nolson /= nolson94) THEN
4292	WRITE(numout,*) nolson,nolson94
4293	CALL ipslerr_p(3,'get_vegcorr', '', '',&
4294	& 'wrong number of OLSON vegetation types.') ! Fatal error
4295	ENDIF !(nolson /= nolson94)
4296
4297	IF (nvm /= nvm13) THEN
4298	WRITE(numout,*) nvm,nvm13
4299	CALL ipslerr_p(3,'get_vegcorr', '', '',&
4300	& 'wrong number of SECHIBA vegetation types.') ! Fatal error
4301	ENDIF !(nvm /= nvm13)
4302
4303	! The carteveg5km cannot be used if the PFTs are not in the standard order
4304	DO ib = 1,nvm
4305	IF (pft_to_mtc(ib) /= ib ) THEN
4306	CALL ipslerr_p(3,'get_vegcorr','You have redefined the order of the 13 PFTS', &
4307	& 'You can not use carteveg5km', 'Use the standard configuration of PFTS' )
4308	ENDIF
4309	ENDDO
4310
4311	!-
4312	! 1 set the indices of non-biospheric surface types to 0.
4313	!-
4314	nobiocorr(:,:) = zero
4315	!-
4316	! 2 Here we construct the correspondance table
4317	! between Olson and the following SECHIBA Classes.
4318	! vegcorr(i,:)+nobiocorr(i,:) = 1. for all i.
4319	!-
4320	! The modified OLSON types found in file carteveg5km.nc
4321	! created by Nicolas Viovy :
4322	! 1 Urban
4323	vegcorr( 1,:) = &
4324	& (/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/)
4325	! 2 Cool low sparse grassland
4326	vegcorr( 2,:) = &
4327	& (/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/)
4328	! 3 Cold conifer forest
4329	vegcorr( 3,:) = &
4330	& (/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/)
4331	! 4 Cold deciduous conifer forest
4332	vegcorr( 4,:) = &
4333	& (/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/)
4334	! 5 Cool Deciduous broadleaf forest
4335	vegcorr( 5,:) = &
4336	& (/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/)
4337	! 6 Cool evergreen broadleaf forests
4338	vegcorr( 6,:) = &
4339	& (/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/)
4340	! 7 Cool tall grasses and shrubs
4341	vegcorr( 7,:) = &
4342	& (/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/)
4343	! 8 Warm C3 tall grasses and shrubs
4344	vegcorr( 8,:) = &
4345	& (/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/)
4346	! 9 Warm C4 tall grases and shrubs
4347	vegcorr( 9,:) = &
4348	& (/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/)
4349	! 10 Bare desert
4350	vegcorr(10,:) = &
4351	& (/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/)
4352	! 11 Cold upland tundra
4353	vegcorr(11,:) = &
4354	& (/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/)
4355	! 12 Cool irrigated grassland
4356	vegcorr(12,:) = &
4357	& (/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/)
4358	! 13 Semi desert
4359	vegcorr(13,:) = &
4360	& (/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/)
4361	! 14 Glacier ice
4362	vegcorr(14,:) = &
4363	& (/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/)
4364	nobiocorr(14,iice) = 1.
4365	! 15 Warm wooded wet swamp
4366	vegcorr(15,:) = &
4367	& (/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/)
4368	! 16 Inland water
4369	vegcorr(16,:) = &
4370	& (/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/)
4371	! 17 sea water
4372	vegcorr(17,:) = &
4373	& (/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/)
4374	! 18 cool shrub evergreen
4375	vegcorr(18,:) = &
4376	& (/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/)
4377	! 19 cold shrub deciduous
4378	vegcorr(19,:) = &
4379	& (/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/)
4380	! 20 Cold evergreen forest and fields
4381	vegcorr(20,:) = &
4382	& (/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/)
4383	! 21 cool rain forest
4384	vegcorr(21,:) = &
4385	& (/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/)
4386	! 22 cold conifer boreal forest
4387	vegcorr(22,:) = &
4388	& (/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/)
4389	! 23 cool conifer forest
4390	vegcorr(23,:) = &
4391	& (/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/)
4392	! 24 warm mixed forest
4393	vegcorr(24,:) = &
4394	& (/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/)
4395	! 25 cool mixed forest
4396	vegcorr(25,:) = &
4397	& (/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/)
4398	! 26 cool broadleaf forest
4399	vegcorr(26,:) = &
4400	& (/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/)
4401	! 27 cool deciduous broadleaf forest
4402	vegcorr(27,:) = &
4403	& (/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/)
4404	! 28 warm montane tropical forest
4405	vegcorr(28,:) = &
4406	& (/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/)
4407	! 29 warm seasonal tropical forest
4408	vegcorr(29,:) = &
4409	& (/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/)
4410	! 30 cool crops and towns
4411	vegcorr(30,:) = &
4412	& (/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/)
4413	! 31 warm crops and towns
4414	vegcorr(31,:) = &
4415	& (/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/)
4416	! 32 cool crops and towns
4417	vegcorr(32,:) = &
4418	& (/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/)
4419	! 33 warm dry tropical woods
4420	vegcorr(33,:) = &
4421	& (/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/)
4422	! 34 warm tropical rain forest
4423	vegcorr(34,:) = &
4424	& (/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/)
4425	! 35 warm tropical degraded forest
4426	vegcorr(35,:) = &
4427	& (/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/)
4428	! 36 warm corn and beans cropland
4429	vegcorr(36,:) = &
4430	& (/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/)
4431	! 37 cool corn and bean cropland
4432	vegcorr(37,:) = &
4433	& (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0/)
4434	! 38 warm rice paddy and field
4435	vegcorr(38,:) = &
4436	& (/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/)
4437	! 39 hot irrigated cropland
4438	vegcorr(39,:) = &
4439	& (/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/)
4440	! 40 cool irrigated cropland
4441	vegcorr(40,:) = &
4442	& (/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/)
4443	! 41 cold irrigated cropland
4444	vegcorr(41,:) = &
4445	& (/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/)
4446	! 42 cool grasses and shrubs
4447	vegcorr(42,:) = &
4448	& (/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/)
4449	! 43 hot and mild grasses and shrubs
4450	vegcorr(43,:) = &
4451	& (/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/)
4452	! 44 cold grassland
4453	vegcorr(44,:) = &
4454	& (/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/)
4455	! 45 Savanna (woods) C3
4456	vegcorr(45,:) = &
4457	& (/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/)
4458	! 46 Savanna woods C4
4459	vegcorr(46,:) = &
4460	& (/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/)
4461	! 47 Mire, bog, fen
4462	vegcorr(47,:) = &
4463	& (/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/)
4464	! 48 Warm marsh wetland
4465	vegcorr(48,:) = &
4466	& (/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/)
4467	! 49 cold marsh wetland
4468	vegcorr(49,:) = &
4469	& (/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/)
4470	! 50 mediteraean scrub
4471	vegcorr(50,:) = &
4472	& (/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/)
4473	! 51 Cool dry woody scrub
4474	vegcorr(51,:) = &
4475	& (/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/)
4476	! 52 Warm dry evergreen woods
4477	vegcorr(52,:) = &
4478	& (/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/)
4479	! 53 Volcanic rocks
4480	vegcorr(53,:) = &
4481	& (/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/)
4482	! 54 sand desert
4483	vegcorr(54,:) = &
4484	& (/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/)
4485	! 55 warm semi desert shrubs
4486	vegcorr(55,:) = &
4487	& (/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/)
4488	! 56 cool semi desert shrubs
4489	vegcorr(56,:) = &
4490	& (/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/)
4491	! 57 semi desert sage
4492	vegcorr(57,:) = &
4493	& (/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/)
4494	! 58 Barren tundra
4495	vegcorr(58,:) = &
4496	& (/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/)
4497	! 59 cool southern hemisphere mixed forest
4498	vegcorr(59,:) = &
4499	& (/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/)
4500	! 60 cool fields and woods
4501	vegcorr(60,:) = &
4502	& (/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/)
4503	! 61 warm forest and filed
4504	vegcorr(61,:) = &
4505	& (/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/)
4506	! 62 cool forest and field
4507	vegcorr(62,:) = &
4508	& (/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/)
4509	! 63 warm C3 fields and woody savanna
4510	vegcorr(63,:) = &
4511	& (/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/)
4512	! 64 warm C4 fields and woody savanna
4513	vegcorr(64,:) = &
4514	& (/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/)
4515	! 65 cool fields and woody savanna
4516	vegcorr(65,:) = &
4517	& (/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/)
4518	! 66 warm succulent and thorn scrub
4519	vegcorr(66,:) = &
4520	& (/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/)
4521	! 67 cold small leaf mixed woods
4522	vegcorr(67,:) = &
4523	& (/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/)
4524	! 68 cold deciduous and mixed boreal fores
4525	vegcorr(68,:) = &
4526	& (/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/)
4527	! 69 cold narrow conifers
4528	vegcorr(69,:) = &
4529	& (/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/)
4530	! 70 cold wooded tundra
4531	vegcorr(70,:) = &
4532	& (/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/)
4533	! 71 cold heath scrub
4534	vegcorr(71,:) = &
4535	& (/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/)
4536	! 72 Polar and alpine desert
4537	vegcorr(72,:) = &
4538	& (/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/)
4539	! 73 warm Mangrove
4540	vegcorr(73,:) = &
4541	& (/0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
4542	! 74 cool crop and water mixtures
4543	vegcorr(74,:) = &
4544	& (/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/)
4545	! 75 cool southern hemisphere mixed forest
4546	vegcorr(75,:) = &
4547	& (/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/)
4548	! 76 cool moist eucalyptus
4549	vegcorr(76,:) = &
4550	& (/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/)
4551	! 77 warm rain green tropical forest
4552	vegcorr(77,:) = &
4553	& (/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/)
4554	! 78 warm C3 woody savanna
4555	vegcorr(78,:) = &
4556	& (/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/)
4557	! 79 warm C4 woody savanna
4558	vegcorr(79,:) = &
4559	& (/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/)
4560	! 80 cool woody savanna
4561	vegcorr(80,:) = &
4562	& (/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/)
4563	! 81 cold woody savanna
4564	vegcorr(81,:) = &
4565	& (/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/)
4566	! 82 warm broadleaf crops
4567	vegcorr(82,:) = &
4568	& (/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/)
4569	! 83 warm C3 grass crops
4570	vegcorr(83,:) = &
4571	& (/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/)
4572	! 84 warm C4 grass crops
4573	vegcorr(84,:) = &
4574	& (/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/)
4575	! 85 cool grass crops
4576	vegcorr(85,:) = &
4577	& (/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/)
4578	! 86 warm C3 crops grass,shrubs
4579	vegcorr(86,:) = &
4580	& (/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/)
4581	! 87 cool crops,grass,shrubs
4582	vegcorr(87,:) = &
4583	& (/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/)
4584	! 88 warm evergreen tree crop
4585	vegcorr(88,:) = &
4586	& (/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/)
4587	! 89 cool evergreen tree crop
4588	vegcorr(89,:) = &
4589	& (/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/)
4590	! 90 cold evergreen tree crop
4591	vegcorr(90,:) = &
4592	& (/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/)
4593	! 91 warm deciduous tree crop
4594	vegcorr(91,:) = &
4595	& (/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/)
4596	! 92 cool deciduous tree crop
4597	vegcorr(92,:) = &
4598	& (/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/)
4599	! 93 cold deciduous tree crop
4600	vegcorr(93,:) = &
4601	& (/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/)
4602	! 94 wet sclerophylic forest
4603	vegcorr(94,:) = &
4604	& (/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/)
4605	!-
4606	! 3 Check the mapping for the Olson types which are going into the
4607	! the veget and nobio array.
4608	!-
4609	DO ib=1,nolson
4610	!
4611	IF ( ABS(SUM(vegcorr(ib,:))+SUM(nobiocorr(ib,:))-1.0) &
4612	& > EPSILON(1.0)) THEN
4613	WRITE(numout,*) 'Wrong correspondance for Olson type :', ib
4614	CALL ipslerr_p(3,'get_vegcorr', '', '',&
4615	& 'Wrong correspondance for Olson type.') ! Fatal error
4616	ENDIF
4617	!
4618	ENDDO ! Loop over the # Olson type
4619
4620
4621	END SUBROUTINE get_vegcorr
4622
4623	!! ================================================================================================================================
4624	!! SUBROUTINE : get_soilcorr_zobler
4625	!!
4626	!>\BRIEF The "get_soilcorr_zobler" routine defines the table of correspondence
4627	!! between the Zobler types and the three texture types known by SECHIBA and STOMATE :
4628	!! silt, sand and clay.
4629	!!
4630	!! DESCRIPTION : get_soilcorr_zobler is needed if you use soils_param.nc .\n
4631	!! The data from this file is then interpolated to the grid of the model. \n
4632	!! The aim is to get fractions for sand loam and clay in each grid box.\n
4633	!! This information is used for soil hydrology and respiration.
4634	!!
4635	!!
4636	!! RECENT CHANGE(S): None
4637	!!
4638	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
4639	!!
4640	!! REFERENCE(S) :
4641	!! - Zobler L., 1986, A World Soil File for global climate modelling. NASA Technical memorandum 87802. NASA
4642	!! Goddard Institute for Space Studies, New York, U.S.A.
4643	!!
4644	!! FLOWCHART : None
4645	!! \n
4646	!_ ================================================================================================================================
4647
4648	SUBROUTINE get_soilcorr_zobler (nzobler,textfrac_table)
4649
4650	IMPLICIT NONE
4651
4652	!! 0. Variables and parameters declaration
4653
4654	INTEGER(i_std),PARAMETER :: nbtypes_zobler = 7 !! Number of Zobler types (unitless)
4655
4656	!! 0.1 Input variables
4657
4658	INTEGER(i_std),INTENT(in) :: nzobler !! Size of the array (unitless)
4659
4660	!! 0.2 Output variables
4661
4662	REAL(r_std),DIMENSION(nzobler,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
4663	!! and granulometric composition (0-1, unitless)
4664
4665	!! 0.4 Local variables
4666
4667	INTEGER(i_std) :: ib !! Indice (unitless)
4668
4669	!_ ================================================================================================================================
4670
4671	!-
4672	! 0. Check consistency
4673	!-
4674	IF (nzobler /= nbtypes_zobler) THEN
4675	CALL ipslerr_p(3,'get_soilcorr_zobler', 'nzobler /= nbtypes_zobler',&
4676	& 'We do not have the correct number of classes', &
4677	& ' in the code for the file.') ! Fatal error
4678	ENDIF
4679
4680	!-
4681	! 1. Textural fraction for : silt sand clay
4682	!-
4683	textfrac_table(1,:) = (/ 0.12, 0.82, 0.06 /)
4684	textfrac_table(2,:) = (/ 0.32, 0.58, 0.10 /)
4685	textfrac_table(3,:) = (/ 0.39, 0.43, 0.18 /)
4686	textfrac_table(4,:) = (/ 0.15, 0.58, 0.27 /)
4687	textfrac_table(5,:) = (/ 0.34, 0.32, 0.34 /)
4688	textfrac_table(6,:) = (/ 0.00, 1.00, 0.00 /)
4689	textfrac_table(7,:) = (/ 0.39, 0.43, 0.18 /)
4690
4691
4692	!-
4693	! 2. Check the mapping for the Zobler types which are going into the ORCHIDEE textures classes
4694	!-
4695	DO ib=1,nzobler ! Loop over # classes soil
4696
4697	IF (ABS(SUM(textfrac_table(ib,:))-1.0) > EPSILON(1.0)) THEN ! The sum of the textural fractions should not exceed 1 !
4698	WRITE(numout,*) &
4699	& 'Error in the correspondence table', &
4700	& ' sum is not equal to 1 in', ib
4701	WRITE(numout,*) textfrac_table(ib,:)
4702	CALL ipslerr_p(3,'get_soilcorr_zobler', 'SUM(textfrac_table(ib,:)) /= 1.0',&
4703	& '', 'Error in the correspondence table') ! Fatal error
4704	ENDIF
4705
4706	ENDDO ! Loop over # classes soil
4707
4708
4709	END SUBROUTINE get_soilcorr_zobler
4710
4711	!! ================================================================================================================================
4712	!! SUBROUTINE : get_soilcorr_usda
4713	!!
4714	!>\BRIEF The "get_soilcorr_usda" routine defines the table of correspondence
4715	!! between the 12 USDA textural classes and their granulometric composition,
4716	!! as % of silt, sand and clay. This is used to further defien clayfraction.
4717	!!
4718	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
4719	!! The data from this file is then interpolated to the grid of the model. \n
4720	!! The aim is to get fractions for sand loam and clay in each grid box.\n
4721	!! This information is used for soil hydrology and respiration.
4722	!! The default map in this case is derived from Reynolds et al 2000, \n
4723	!! at the 1/12deg resolution, with indices that are consistent with the \n
4724	!! textures tabulated below
4725	!!
4726	!! RECENT CHANGE(S): Created by A. Ducharne on July 02, 2014
4727	!!
4728	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
4729	!!
4730	!! REFERENCE(S) :
4731	!!
4732	!! FLOWCHART : None
4733	!! \n
4734	!_ ================================================================================================================================
4735
4736	SUBROUTINE get_soilcorr_usda (nusda,textfrac_table)
4737
4738	IMPLICIT NONE
4739
4740	!! 0. Variables and parameters declaration
4741
4742	!! 0.1 Input variables
4743
4744	INTEGER(i_std),INTENT(in) :: nusda !! Size of the array (unitless)
4745
4746	!! 0.2 Output variables
4747
4748	REAL(r_std),DIMENSION(nusda,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
4749	!! and granulometric composition (0-1, unitless)
4750
4751	!! 0.4 Local variables
4752
4753	INTEGER(i_std),PARAMETER :: nbtypes_usda = 12 !! Number of USDA texture classes (unitless)
4754	INTEGER(i_std) :: n !! Index (unitless)
4755
4756	!_ ================================================================================================================================
4757
4758	!-
4759	! 0. Check consistency
4760	!-
4761	IF (nusda /= nbtypes_usda) THEN
4762	CALL ipslerr_p(3,'get_soilcorr_usda', 'nusda /= nbtypes_usda',&
4763	& 'We do not have the correct number of classes', &
4764	& ' in the code for the file.') ! Fatal error
4765	ENDIF
4766
4767	!! Parameters for soil type distribution :
4768	!! Sand, Loamy Sand, Sandy Loam, Silt Loam, Silt, Loam, Sandy Clay Loam, Silty Clay Loam, Clay Loam, Sandy Clay, Silty Clay, Clay
4769	! The order comes from constantes_soil.f90
4770	! The corresponding granulometric composition comes from Carsel & Parrish, 1988
4771
4772	!-
4773	! 1. Textural fractions for : sand, clay
4774	!-
4775	textfrac_table(1,2:3) = (/ 0.93, 0.03 /) ! Sand
4776	textfrac_table(2,2:3) = (/ 0.81, 0.06 /) ! Loamy Sand
4777	textfrac_table(3,2:3) = (/ 0.63, 0.11 /) ! Sandy Loam
4778	textfrac_table(4,2:3) = (/ 0.17, 0.19 /) ! Silt Loam
4779	textfrac_table(5,2:3) = (/ 0.06, 0.10 /) ! Silt
4780	textfrac_table(6,2:3) = (/ 0.40, 0.20 /) ! Loam
4781	textfrac_table(7,2:3) = (/ 0.54, 0.27 /) ! Sandy Clay Loam
4782	textfrac_table(8,2:3) = (/ 0.08, 0.33 /) ! Silty Clay Loam
4783	textfrac_table(9,2:3) = (/ 0.30, 0.33 /) ! Clay Loam
4784	textfrac_table(10,2:3) = (/ 0.48, 0.41 /) ! Sandy Clay
4785	textfrac_table(11,2:3) = (/ 0.06, 0.46 /) ! Silty Clay
4786	textfrac_table(12,2:3) = (/ 0.15, 0.55 /) ! Clay
4787
4788	! Fraction of silt
4789
4790	DO n=1,nusda
4791	textfrac_table(n,1) = 1. - textfrac_table(n,2) - textfrac_table(n,3)
4792	END DO
4793
4794	END SUBROUTINE get_soilcorr_usda
4795
4796	!! ================================================================================================================================
4797	!! FUNCTION : tempfunc
4798	!!
4799	!>\BRIEF ! This function interpolates value between ztempmin and ztempmax
4800	!! used for lai detection.
4801	!!
4802	!! DESCRIPTION : This subroutine calculates a scalar between 0 and 1 with the following equation :\n
4803	!! \latexonly
4804	!! \input{constantes_veg_tempfunc.tex}
4805	!! \endlatexonly
4806	!!
4807	!! RECENT CHANGE(S): None
4808	!!
4809	!! RETURN VALUE : tempfunc_result
4810	!!
4811	!! REFERENCE(S) : None
4812	!!
4813	!! FLOWCHART : None
4814	!! \n
4815	!_ ================================================================================================================================
4816
4817	FUNCTION tempfunc (temp_in) RESULT (tempfunc_result)
4818
4819
4820	!! 0. Variables and parameters declaration
4821
4822	REAL(r_std),PARAMETER :: ztempmin=273._r_std !! Temperature for laimin (K)
4823	REAL(r_std),PARAMETER :: ztempmax=293._r_std !! Temperature for laimax (K)
4824	REAL(r_std) :: zfacteur !! Interpolation factor (K^{-2})
4825
4826	!! 0.1 Input variables
4827
4828	REAL(r_std),INTENT(in) :: temp_in !! Temperature (K)
4829
4830	!! 0.2 Result
4831
4832	REAL(r_std) :: tempfunc_result !! (unitless)
4833
4834	!_ ================================================================================================================================
4835
4836	!! 1. Define a coefficient
4837	zfacteur = un/(ztempmax-ztempmin)**2
4838
4839	!! 2. Computes tempfunc
4840	IF (temp_in > ztempmax) THEN
4841	tempfunc_result = un
4842	ELSEIF (temp_in < ztempmin) THEN
4843	tempfunc_result = zero
4844	ELSE
4845	tempfunc_result = un-zfacteur(ztempmax-temp_in)*2
4846	ENDIF !(temp_in > ztempmax)
4847
4848
4849	END FUNCTION tempfunc
4850
4851	!!
4852	!================================================================================================================================
4853	!! SUBROUTINE : slowproc_read_veg_restfile
4854	!!
4855	!>\BRIEF read a vegetation map from a restart file, and interpolate if
4856	!necessary to current model grid
4857	!!
4858	!! DESCRIPTION : (definitions, functional, design, flags):
4859	!!
4860	!! RECENT CHANGE(S): None
4861	!!
4862	!! MAIN OUTPUT VARIABLE(S): :: none
4863	!!
4864	!! REFERENCE(S) : None
4865	!!
4866	!! FLOWCHART : None
4867	!! \n
4868	!_
4869	!================================================================================================================================
4870
4871	SUBROUTINE slowproc_read_veg_restfile(nbpt, lalo, neighbours, resolution, contfrac, veget, frac_nobio )
4872	!
4873	!
4874	!
4875	! 0.1 INPUT
4876	!
4877	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
4878	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
4879	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
4880	! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
4881	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
4882	REAL(r_std),DIMENSION (nbpt), INTENT (in) :: contfrac !! Fraction of continent in the grid
4883	!
4884	! 0.2 OUTPUT
4885	!
4886	REAL(r_std), INTENT(out) :: veget(nbpt,nvm) ! Vegetation fractions
4887	REAL(r_std), INTENT(out) :: frac_nobio(nbpt,nnobio) ! Fraction of the mesh which is covered by ice, lakes, ...
4888	!
4889	! 0.3 LOCAL
4890	!
4891	!
4892	CHARACTER(LEN=80) :: filename
4893	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, jp, it, jj, jv, im, jm
4894	INTEGER(i_std) :: il, ils, ip, ix, iy, imin, jmin, ier
4895	REAL(r_std) :: dlon, dlonmin, dlat, dlatmin
4896	REAL(r_std), PARAMETER :: maxmargin = 5.
4897	!
4898	! 0.4 allocatable
4899	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: veget_max_in
4900	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: frac_nobio_in
4901	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: xx, yy
4902	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: x, y
4903
4904	! get filename
4905	CALL getin_p('vegetation_map_restartfile_name', filename)
4906
4907	! load the sechiba restart file, and read in values of veget_max and
4908	! frac_nobio
4909	CALL flininfo(filename, iml, jml, lml, tml, fid)
4910	!! for now, assume hard-coded 96 x 72 gridpoints because iml and jml
4911	!readings are wrong
4912	im = 96
4913	jm = 72
4914	ALLOCATE(xx(im, jm))
4915	ALLOCATE(yy(im,jm))
4916	ALLOCATE (x(im),y(jm))
4917	ALLOCATE(veget_max_in(im,jm,nvm))
4918	ALLOCATE(frac_nobio_in(im,jm))
4919	WRITE(,) 'cdk debug: '
4920	WRITE(,) 'filename: ', filename
4921	WRITE(,) 'im, jm, lml, tml, fid: ', im, jm, lml, tml, fid
4922
4923	CALL flinget(fid, 'nav_lon', im, jm, 1, 1, 1, 1, xx)
4924	WRITE(,) 'cdk debug read nav_lon '
4925	CALL flinget(fid, 'nav_lat', im, jm, 1, 1, 1, 1, yy)
4926	WRITE(,) 'cdk debug read nav_lat '
4927	CALL flinget(fid, 'veget_max', im, jm, nvm, 1, 1, 1, veget_max_in) !! these need to be double-precision
4928	WRITE(,) 'cdk debug read veget_max '
4929	CALL flinget(fid, 'frac_nobio', im, jm, 1, 1, 1, 1, frac_nobio_in) !! these need to be double-precision
4930	WRITE(,) 'cdk debug read frac_nobio '
4931
4932	x(:) = xx(:,1)
4933	y(:) = yy(1,:)
4934
4935	! prendre la valeur la plus proche
4936	DO ip = 1, nbpt
4937	dlonmin = HUGE(1.)
4938	DO ix = 1,im
4939	dlon = MIN( ABS(lalo(ip,2)-x(ix)), ABS(lalo(ip,2)+360.-x(ix)), ABS(lalo(ip,2)-360.-x(ix)) )
4940	IF ( dlon .LT. dlonmin ) THEN
4941	imin = ix
4942	dlonmin = dlon
4943	ENDIF
4944	ENDDO
4945	dlatmin = HUGE(1.)
4946	DO iy = 1,jm
4947	dlat = ABS(lalo(ip,1)-y(iy))
4948	IF ( dlat .LT. dlatmin ) THEN
4949	jmin = iy
4950	dlatmin = dlat
4951	ENDIF
4952	ENDDO
4953	WRITE(,) 'cdk debug found: ip, dlonmin, dlatmin: ', ip, dlonmin,dlatmin
4954	! imin, jmin, x(imin), y(jmin):,
4955	! if nothing is close, then set as all nobio
4956	IF ( ( dlonmin .LE. maxmargin ) .AND. ( dlatmin .LE. maxmargin ) ) THEN
4957	veget(ip,:) = veget_max_in(imin,jmin,:)
4958	frac_nobio(ip,:) = frac_nobio_in(imin,jmin)
4959	ELSE
4960	veget(ip,:) = zero
4961	frac_nobio(ip,:) = un
4962	ENDIF
4963	WRITE(,) 'cdk debug: ip, veget, frac_nobio: ', veget(ip,:),frac_nobio(ip,:)
4964
4965	ENDDO
4966
4967	DEALLOCATE(xx)
4968	DEALLOCATE(yy)
4969	DEALLOCATE(x)
4970	DEALLOCATE(y)
4971	DEALLOCATE(veget_max_in)
4972	DEALLOCATE(frac_nobio_in)
4973	!
4974	RETURN
4975	!
4976	END SUBROUTINE slowproc_read_veg_restfile
4977
4978
4979	!! ================================================================================================================================
4980	!! SUBROUTINE : slowproc_checkveget
4981	!!
4982	!>\BRIEF To verify the consistency of the various fractions defined within the grid box after having been
4983	!! been updated by STOMATE or the standard procedures.
4984	!!
4985	!! DESCRIPTION : (definitions, functional, design, flags):
4986	!!
4987	!! RECENT CHANGE(S): None
4988	!!
4989	!! MAIN OUTPUT VARIABLE(S): :: none
4990	!!
4991	!! REFERENCE(S) : None
4992	!!
4993	!! FLOWCHART : None
4994	!! \n
4995	!_ ================================================================================================================================
4996	!
4997	SUBROUTINE slowproc_checkveget(nbpt, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
4998
4999	! 0.1 INPUT
5000	!
5001	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
5002	REAL(r_std),DIMENSION (nbpt,nnobio), INTENT(in) :: frac_nobio ! Fraction of ice,lakes,cities, ... (unitless)
5003	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget_max ! Maximum fraction of vegetation type including none biological fraction (unitless)
5004	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget ! Vegetation fractions
5005	REAL(r_std),DIMENSION (nbpt), INTENT(in) :: tot_bare_soil ! Total evaporating bare soil fraction within the mesh
5006	REAL(r_std),DIMENSION (nbpt,nstm), INTENT(in) :: soiltile ! Fraction of soil tiles in the gridbox (unitless)
5007
5008	! 0.3 LOCAL
5009	!
5010	INTEGER(i_std) :: ji, jn, jv
5011	REAL(r_std) :: epsilocal !! A very small value
5012	REAL(r_std) :: totfrac
5013	CHARACTER(len=80) :: str1, str2
5014
5015	!_ ================================================================================================================================
5016
5017	!
5018	! There is some margin added as the computing errors might bring us above EPSILON(un)
5019	!
5020	epsilocal = EPSILON(un)*1000.
5021
5022	!! 1.0 Verify that none of the fractions are smaller than min_vegfrac, without beeing zero.
5023	!!
5024	DO ji=1,nbpt
5025	DO jn=1,nnobio
5026	IF ( frac_nobio(ji,jn) > epsilocal .AND. frac_nobio(ji,jn) < min_vegfrac ) THEN
5027	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
5028	WRITE(str2,'("The small value obtained is ", E14.4)') frac_nobio(ji,jn)
5029	CALL ipslerr_p (3,'slowproc_checkveget', &
5030	"frac_nobio is larger than zero but smaller than min_vegfrac.", str1, str2)
5031	ENDIF
5032	ENDDO
5033	END DO
5034
5035	IF (.NOT. ok_dgvm) THEN
5036	DO ji=1,nbpt
5037	DO jv=1,nvm
5038	IF ( veget_max(ji,jv) > epsilocal .AND. veget_max(ji,jv) < min_vegfrac ) THEN
5039	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
5040	WRITE(str2,'("The small value obtained is ", E14.4)') veget_max(ji,jv)
5041	CALL ipslerr_p (3,'slowproc_checkveget', &
5042	"veget_max is larger than zero but smaller than min_vegfrac.", str1, str2)
5043	ENDIF
5044	ENDDO
5045	ENDDO
5046	END IF
5047
5048	!! 2.0 verify that with all the fractions we cover the entire grid box
5049	!!
5050	DO ji=1,nbpt
5051	totfrac = zero
5052	DO jn=1,nnobio
5053	totfrac = totfrac + frac_nobio(ji,jn)
5054	ENDDO
5055	DO jv=1,nvm
5056	totfrac = totfrac + veget_max(ji,jv)
5057	ENDDO
5058	IF ( ABS(totfrac - un) > epsilocal) THEN
5059	WRITE(str1,'("This occurs on grid box", I8)') ji
5060	WRITE(str2,'("The sum over all fraction and error are ", E14.4, E14.4)') totfrac, ABS(totfrac - un)
5061	CALL ipslerr_p (3,'slowproc_checkveget', &
5062	"veget_max + frac_nobio is not equal to 1.", str1, str2)
5063	WRITE(,) "EPSILON =", epsilocal
5064	ENDIF
5065	ENDDO
5066
5067	!! 3.0 Verify that veget is smaller or equal to veget_max
5068	!!
5069	DO ji=1,nbpt
5070	DO jv=1,nvm
5071	IF ( jv == ibare_sechiba ) THEN
5072	IF ( ABS(veget(ji,jv) - veget_max(ji,jv)) > epsilocal ) THEN
5073	WRITE(str1,'("This occurs on grid box", I8)') ji
5074	WRITE(str2,'("The difference is ", E14.4)') veget(ji,jv) - veget_max(ji,jv)
5075	CALL ipslerr_p (3,'slowproc_checkveget', &
5076	"veget is not equal to veget_max on bare soil.", str1, str2)
5077	ENDIF
5078	ELSE
5079	IF ( veget(ji,jv) > veget_max(ji,jv) ) THEN
5080	WRITE(str1,'("This occurs on grid box", I8)') ji
5081	WRITE(str2,'("The values for veget and veget_max :", F8.4, F8.4)') veget(ji,jv), veget_max(ji,jv)
5082	CALL ipslerr_p (3,'slowproc_checkveget', &
5083	"veget is greater than veget_max.", str1, str2)
5084	ENDIF
5085	ENDIF
5086	ENDDO
5087	ENDDO
5088
5089	!! 4.0 Test tot_bare_soil in relation to the other variables
5090	!!
5091	DO ji=1,nbpt
5092	totfrac = zero
5093	DO jv=1,nvm
5094	totfrac = totfrac + (veget_max(ji,jv) - veget(ji,jv))
5095	ENDDO
5096	! add the bare soil fraction to totfrac
5097	totfrac = totfrac + veget(ji,ibare_sechiba)
5098	! do the test
5099	IF ( ABS(totfrac - tot_bare_soil(ji)) > epsilocal ) THEN
5100	WRITE(str1,'("This occurs on grid box", I8)') ji
5101	WRITE(str2,'("The values for tot_bare_soil, tot frac and error :", F8.4, F8.4, E14.4)') &
5102	& tot_bare_soil(ji), totfrac, ABS(totfrac - tot_bare_soil(ji))
5103	CALL ipslerr_p (3,'slowproc_checkveget', &
5104	"tot_bare_soil does not correspond to the total bare soil fraction.", str1, str2)
5105	ENDIF
5106	ENDDO
5107
5108	!! 5.0 Test that soiltile has the right sum
5109	!!
5110	DO ji=1,nbpt
5111	totfrac = SUM(soiltile(ji,:))
5112	IF ( ABS(totfrac - un) > epsilocal ) THEN
5113	WRITE(numout,*) "soiltile does not sum-up to one. This occurs on grid box", ji
5114	WRITE(numout,*) "The soiltile for ji are :", soiltile(ji,:)
5115	CALL ipslerr_p (2,'slowproc_checkveget', &
5116	"soiltile does not sum-up to one.", "", "")
5117	ENDIF
5118	ENDDO
5119
5120	END SUBROUTINE slowproc_checkveget
5121
5122	!! ================================================================================================================================
5123	!! SUBROUTINE : slowproc_change_frac
5124	!!
5125	!>\BRIEF Update the vegetation fractions
5126	!!
5127	!! DESCRIPTION : Update the vegetation fractions. This subroutine is called in the same time step as lcchange in stomatelpj has
5128	!! has been done. This subroutine is called after the diagnostics have been written in sechiba_main.
5129	!!
5130	!! RECENT CHANGE(S): None
5131	!!
5132	!! MAIN OUTPUT VARIABLE(S): :: veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile
5133	!!
5134	!! REFERENCE(S) : None
5135	!!
5136	!! FLOWCHART : None
5137	!! \n
5138	!_ ================================================================================================================================
5139
5140	SUBROUTINE slowproc_change_frac(kjpindex, f_rot_sech, lai, &
5141	veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile)
5142	!
5143	! 0. Declarations
5144	!
5145	! 0.1 Input variables
5146	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
5147	LOGICAL,DIMENSION(kjpindex),INTENT(in) :: f_rot_sech !! whether a grid is under rotation
5148	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m^2 m^{-2})
5149
5150	! 0.2 Output variables
5151	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
5152	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget !! Fraction of vegetation type in the mesh (unitless)
5153	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
5154	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
5155	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
5156	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
5157
5158	! 0.3 Local variables
5159	INTEGER(i_std) :: ji, jv !! Loop index
5160
5161
5162	!! Update vegetation fractions with the values coming from the vegetation file read in slowproc_readvegetmax.
5163	!! Partial update has been taken into account for the case with DGVM and AGRICULTURE in slowproc_readvegetmax.
5164	veget_max = veget_max_new
5165	frac_nobio = frac_nobio_new
5166
5167	!! Verification and correction on veget_max, calculation of veget and soiltile.
5168	CALL slowproc_veget (kjpindex, f_rot_sech, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile)
5169
5170	!! Calculate tot_bare_soil needed in hydrol, diffuco and condveg (fraction of bare soil in the mesh)
5171	tot_bare_soil(:) = veget_max(:,1)
5172	DO jv = 2, nvm
5173	DO ji =1, kjpindex
5174	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
5175	ENDDO
5176	END DO
5177
5178	!! Do some basic tests on the surface fractions updated above
5179	CALL slowproc_checkveget(kjpindex, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
5180
5181	END SUBROUTINE slowproc_change_frac
5182
5183	!spitfire
5184	SUBROUTINE slowproc_read_data(nbpt, lalo, resolution, proxydata,data_filename,field_name)
5185
5186	!
5187	!
5188	!
5189	! 0.1 INPUT
5190	!
5191	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
5192	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
5193	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
5194	!
5195	! 0.2 OUTPUT
5196	!
5197	REAL(r_std), INTENT(out) :: proxydata(:,:) ! lightn read variable and re-dimensioned
5198	!
5199	! 0.3 LOCAL
5200	!
5201	REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8
5202	!
5203	!
5204	CHARACTER(LEN=*),INTENT(in) :: data_filename
5205	CHARACTER(LEN=*),INTENT(in) :: field_name
5206	INTEGER(i_std) :: iml, jml, ijml, i, j, ik, lml, tml, fid, ib, jb,ip, jp, vid, ai,iki,jkj
5207	INTEGER(i_std) :: nb_coord,nb_dim,nb_var,nb_gat
5208	LOGICAL :: l_ex
5209	INTEGER,DIMENSION(1) :: l_d_w
5210	REAL(r_std) :: lev(1), date, dt, coslat, pi
5211	INTEGER(i_std) :: itau(1)
5212
5213	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu
5214	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu, mask
5215	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful
5216	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lightn_orig
5217	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: lightn_lu !LOOP: as 4D needed?
5218	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
5219	INTEGER, DIMENSION(nbpt) :: n_origlightn
5220	INTEGER, DIMENSION(nbpt) :: n_found
5221
5222	CHARACTER(LEN=80) :: meter
5223	REAL(r_std) :: prog, sumf
5224	LOGICAL :: found
5225	INTEGER :: idi,jdi, ilast, jlast, jj, ii, jv, inear, iprog
5226	REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
5227	!
5228	pi = 4. * ATAN(1.)
5229	!
5230	!Config Key = Lightn_FILE
5231	!Config Desc = Name of file from which the lightn climatology is to be read
5232	!Config If = !lightn
5233	!Config Def = lightn_climatology_otd.nc
5234	!Config Help = The name of the file to be opened to read the lightnin flash rate
5235	!Config map is to be given here. Resolution is 1.0x1.0deg.
5236
5237	!
5238	CALL flininfo(data_filename, iml, jml, lml, tml, fid)
5239
5240	!
5241	!
5242	ALLOCATE(lon_lu(iml))
5243	ALLOCATE(lat_lu(jml))
5244	ALLOCATE(lightn_lu(iml,jml,tml))
5245	ALLOCATE(mask_lu(iml,jml))
5246	!
5247	!
5248	WRITE(numout,*) 'input filename : ', data_filename
5249	CALL flinget(fid, 'lon', iml, 0, 0, 0, 1, 1, lon_lu)
5250	CALL flinget(fid, 'lat', jml, 0, 0, 0, 1, 1, lat_lu)
5251	CALL flinget(fid, field_name, iml, jml, 0, tml, 1, tml, lightn_lu)
5252	!
5253	WRITE(numout,*) 'cordinate information: ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml),tml
5254
5255	!
5256	!
5257	ijml=iml*jml
5258	ALLOCATE(lon_ful(ijml))
5259	ALLOCATE(lat_ful(ijml))
5260	ALLOCATE(lightn_orig(ijml,tml))
5261	ALLOCATE(mask(ijml))
5262
5263
5264	DO i=1,iml
5265	DO j=1,jml
5266	iki=(j-1)*iml+i
5267	lon_ful(iki)=lon_lu(i)
5268	lat_ful(iki)=lat_lu(j)
5269	lightn_orig(iki,:)=lightn_lu(i,j,:)
5270	IF (lightn_lu(i,j,1).gt.-9000.) &
5271	mask(iki) = 1.0
5272	ENDDO
5273	ENDDO
5274
5275	!
5276	ALLOCATE(lon_up(nbpt))
5277	ALLOCATE(lon_low(nbpt))
5278	ALLOCATE(lat_up(nbpt))
5279	ALLOCATE(lat_low(nbpt))
5280	!
5281	DO ib =1, nbpt
5282	!
5283	! We find the 4 limits of the grid-box. As we transform the resolution of the model
5284	! into longitudes and latitudes we do not have the problem of periodicity.
5285	! coslat is a help variable here !
5286	!
5287	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
5288	!
5289	lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
5290	lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
5291	!
5292	coslat = pi/180. * R_Earth
5293	!
5294	lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
5295	lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
5296	!
5297	!
5298	!
5299	ENDDO
5300	!
5301	! Get the limits of the integration domaine so that we can speed up the calculations
5302	!
5303	domaine_lon_min = MINVAL(lon_low)
5304	domaine_lon_max = MAXVAL(lon_up)
5305	domaine_lat_min = MINVAL(lat_low)
5306	domaine_lat_max = MAXVAL(lat_up)
5307	!
5308	! Ensure that the fine grid covers the whole domain
5309	WHERE ( lon_ful(:) .LT. domaine_lon_min )
5310	lon_ful(:) = lon_ful(:) + 360.
5311	ENDWHERE
5312	!
5313	WHERE ( lon_ful(:) .GT. domaine_lon_max )
5314	lon_ful(:) = lon_ful(:) - 360.
5315	ENDWHERE
5316	!
5317	WRITE(numout,*) 'Interpolating the --',field_name,'-- map :'
5318	WRITE(numout,'(2a40)')'0%--------------------------------------', &
5319	& '------------------------------------100%'
5320	!
5321	ilast = 1
5322	n_origlightn(:) = 0.
5323	proxydata(:,:) = 0.
5324	!
5325	DO ip=1,ijml
5326	!
5327	! Give a progress meter
5328	!
5329	iprog = NINT(float(ip)/float(ijml)79.) - NINT(float(ip-1)/float(ijml)79.)
5330	IF ( iprog .NE. 0 ) THEN
5331	WRITE(numout,'(a1,$)') 'y'
5332	ENDIF
5333	!
5334	! Only start looking for its place in the smaler grid if we are within the domaine
5335	! That should speed up things !
5336	!
5337	IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
5338	( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
5339	( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
5340	( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
5341	!
5342	! look for point on GCM grid which this point on fine grid belongs to.
5343	! First look at the point on the model grid where we arrived just before. There is
5344	! a good chance that neighbouring points on the fine grid fall into the same model
5345	! grid box.
5346	!
5347	IF ( ( lon_ful(ip) .GE. lon_low(ilast) ) .AND. &
5348	( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
5349	( lat_ful(ip) .GE. lat_low(ilast) ) .AND. &
5350	( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
5351	!
5352	! We were lucky
5353	!
5354	IF (mask(ip) .GT. 0) THEN
5355	n_origlightn(ilast) = n_origlightn(ilast) + 1
5356	DO j=1,tml
5357	proxydata(ilast,j) = proxydata(ilast,j) + lightn_orig(ip,j)
5358	ENDDO
5359	ENDIF
5360	!
5361	ELSE
5362	!
5363	! Otherwise, look everywhere.
5364	! Begin close to last grid point.
5365	!
5366	found = .FALSE.
5367	idi = 1
5368	!
5369	DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
5370
5371	!
5372	! forward and backward
5373	!
5374	DO ii = 1,2
5375	!
5376	IF ( ii .EQ. 1 ) THEN
5377	ib = ilast - idi
5378	ELSE
5379	ib = ilast + idi
5380	ENDIF
5381	!
5382	IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
5383	IF ( ( lon_ful(ip) .GE. lon_low(ib) ) .AND. &
5384	( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
5385	( lat_ful(ip) .GE. lat_low(ib) ) .AND. &
5386	( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
5387	!
5388	IF (mask(ip) .gt. 0) THEN
5389	! DO i=1,nvm
5390	DO j=1,tml
5391	proxydata(ib,j) = proxydata(ib,j) + lightn_orig(ip,j)
5392	ENDDO
5393	! ENDDO
5394	n_origlightn(ib) = n_origlightn(ib) + 1
5395	ENDIF
5396	ilast = ib
5397	found = .TRUE.
5398	!
5399	ENDIF
5400	ENDIF
5401	!
5402	ENDDO
5403	!
5404	idi = idi + 1
5405	!
5406	ENDDO
5407	!
5408	ENDIF ! lucky/not lucky
5409	!
5410	ENDIF ! in the domain
5411	ENDDO
5412
5413	WRITE(numout,*) ''
5414
5415	! determine fraction of lightning points in each box of the coarse grid
5416	DO ip=1,nbpt
5417	IF ( n_origlightn(ip) .GT. 0 ) THEN
5418	proxydata(ip,:) = proxydata(ip,:)/REAL(n_origlightn(ip),r_std)
5419	ELSE
5420	!
5421	WRITE(numout,*) 'PROBLEM, no point in the ba map found for this grid box'
5422	WRITE(numout,*) 'Longitude range : ', lon_low(ip), lon_up(ip)
5423	WRITE(numout,*) 'Latitude range : ', lat_low(ip), lat_up(ip)
5424	!
5425	WRITE(numout,*) 'Looking for nearest point on the ba map file'
5426	CALL slowproc_nearest (ijml, lon_ful, lat_ful, &
5427	lalo(ip,2), lalo(ip,1), inear)
5428	WRITE(numout,*) 'Coordinates of the nearest point, ',inear,' :', &
5429	lon_ful(inear),lat_ful(inear)
5430	proxydata(ip,:) = lightn_orig(inear,:)
5431	ENDIF
5432	ENDDO
5433	!
5434	IF (printlev_loc>=1) WRITE(numout,*) '--',field_name,'--: Interpolation Done'
5435	WRITE(numout,*) ''
5436
5437	!
5438	DEALLOCATE(lon_up)
5439	DEALLOCATE(lon_low)
5440	DEALLOCATE(lat_up)
5441	DEALLOCATE(lat_low)
5442	DEALLOCATE(lat_ful)
5443	DEALLOCATE(lon_ful)
5444	DEALLOCATE(lat_lu)
5445	DEALLOCATE(lon_lu)
5446	DEALLOCATE(lightn_lu)
5447	DEALLOCATE(lightn_orig)
5448	DEALLOCATE(mask)
5449	!
5450	RETURN
5451	!
5452	END SUBROUTINE slowproc_read_data
5453
5454	!--LOOP
5455	SUBROUTINE slowproc_read_annual(nbpt, lalo, resolution, popd,data_filename,field_name)
5456
5457	!
5458	!
5459	!
5460	! 0.1 INPUT
5461	!
5462	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
5463	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
5464	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
5465	!
5466	! 0.2 OUTPUT
5467	!
5468	REAL(r_std), INTENT(out) :: popd(nbpt) ! population density read variable and re-dimensioned
5469	!
5470	! 0.3 LOCAL
5471	!
5472	REAL(r_std), PARAMETER :: R_Earth = 6378000., min_sechiba=1.E-8
5473	!
5474	!
5475	CHARACTER(LEN=*),INTENT(in) :: data_filename
5476	CHARACTER(LEN=*),INTENT(in) :: field_name
5477	INTEGER(i_std) :: iml, jml, ijml, i, j, ik, lml, tml, fid, ib, jb,ip, jp, vid, ai,iki,jkj
5478	INTEGER(i_std) :: nb_coord,nb_dim,nb_var,nb_gat
5479	LOGICAL :: l_ex
5480	INTEGER,DIMENSION(1) :: l_d_w
5481	REAL(r_std) :: lev(1), date, dt, coslat, pi
5482	INTEGER(i_std) :: itau(1)
5483
5484	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_lu
5485	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu, mask
5486	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_ful, lon_ful
5487	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: popden_orig
5488	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: popden_lu !LOOP: as 4D needed?
5489	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_up, lon_low, lat_up, lat_low
5490	INTEGER, DIMENSION(nbpt) :: n_origpopden
5491	INTEGER, DIMENSION(nbpt) :: n_found
5492
5493	CHARACTER(LEN=80) :: meter
5494	REAL(r_std) :: prog, sumf
5495	LOGICAL :: found
5496	INTEGER :: idi,jdi, ilast, jlast, jj, ii, jv, inear, iprog
5497	REAL(r_std) :: domaine_lon_min, domaine_lon_max, domaine_lat_min, domaine_lat_max
5498	!
5499	pi = 4. * ATAN(1.)
5500
5501	!
5502	!Config Key = POPDENS_FILE
5503	!Config Desc = Name of file from which the popden climatology is to be read
5504	!Config If = !popden
5505	!Config Def = popdens.nc
5506	!Config Help = The name of the file to be opened to read the population density
5507	!Config map is to be given here. Resolution is 1.0x1.0deg.
5508
5509	!
5510	CALL flininfo(data_filename, iml, jml, lml, tml, fid)
5511
5512	!
5513	!
5514	ALLOCATE(lon_lu(iml))
5515	ALLOCATE(lat_lu(jml))
5516	ALLOCATE(popden_lu(iml,jml))
5517	ALLOCATE(mask_lu(iml,jml))
5518	!
5519	!
5520	WRITE(numout,*) 'input filename : ', data_filename
5521	CALL flinget(fid, 'lon', iml, 0, 0, 0, 1, 1, lon_lu)
5522	CALL flinget(fid, 'lat', jml, 0, 0, 0, 1, 1, lat_lu)
5523	CALL flinget(fid, field_name, iml, jml, 0, tml, 1, 1, popden_lu)
5524
5525	!
5526	WRITE(numout,*) 'cordinate information: ', lon_lu(1), lon_lu(iml),lat_lu(1), lat_lu(jml),tml
5527
5528	!
5529	!
5530	ijml=iml*jml
5531	ALLOCATE(lon_ful(ijml))
5532	ALLOCATE(lat_ful(ijml))
5533	ALLOCATE(popden_orig(ijml))
5534	ALLOCATE(mask(ijml))
5535
5536
5537	DO i=1,iml
5538	DO j=1,jml
5539	iki=(j-1)*iml+i
5540	lon_ful(iki)=lon_lu(i)
5541	lat_ful(iki)=lat_lu(j)
5542	popden_orig(iki)=popden_lu(i,j)
5543	IF (popden_lu(i,j).gt.-9000.) &
5544	mask(iki) = 1.0
5545	ENDDO
5546	ENDDO
5547
5548	!
5549	WHERE ( popden_orig(:) .LT. 0 )
5550	popden_orig(:) = 0.
5551	ENDWHERE
5552	!
5553	!
5554	ALLOCATE(lon_up(nbpt))
5555	ALLOCATE(lon_low(nbpt))
5556	ALLOCATE(lat_up(nbpt))
5557	ALLOCATE(lat_low(nbpt))
5558	!
5559	DO ib =1, nbpt
5560	!
5561	! We find the 4 limits of the grid-box. As we transform the resolution of the model
5562	! into longitudes and latitudes we do not have the problem of periodicity.
5563	! coslat is a help variable here !
5564	!
5565	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
5566	!
5567	lon_up(ib) = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
5568	lon_low(ib) = lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
5569	!
5570	coslat = pi/180. * R_Earth
5571	!
5572	lat_up(ib) = lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
5573	lat_low(ib) = lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
5574	!
5575	!
5576	!
5577	ENDDO
5578	!
5579	! Get the limits of the integration domaine so that we can speed up the calculations
5580	!
5581	domaine_lon_min = MINVAL(lon_low)
5582	domaine_lon_max = MAXVAL(lon_up)
5583	domaine_lat_min = MINVAL(lat_low)
5584	domaine_lat_max = MAXVAL(lat_up)
5585	!
5586	!
5587	! Ensure that the fine grid covers the whole domain
5588	WHERE ( lon_ful(:) .LT. domaine_lon_min )
5589	lon_ful(:) = lon_ful(:) + 360.
5590	ENDWHERE
5591	!
5592	WHERE ( lon_ful(:) .GT. domaine_lon_max )
5593	lon_ful(:) = lon_ful(:) - 360.
5594	ENDWHERE
5595	!
5596	WRITE(numout,*) 'Interpolating the --',field_name,'-- map :'
5597	WRITE(numout,'(2a40)')'0%--------------------------------------', &
5598	& '------------------------------------100%'
5599	!
5600	ilast = 1
5601	n_origpopden(:) = 0
5602	popd(:) = 0.
5603	!
5604	DO ip=1,ijml
5605	!
5606	! Give a progress meter
5607	!
5608	iprog = NINT(float(ip)/float(ijml)79.) - NINT(float(ip-1)/float(ijml)79.)
5609	IF ( iprog .NE. 0 ) THEN
5610	WRITE(numout,'(a1,$)') 'y'
5611	ENDIF
5612	!
5613	! Only start looking for its place in the smaler grid if we are within the domaine
5614	! That should speed up things !
5615	!
5616	IF ( ( lon_ful(ip) .GE. domaine_lon_min ) .AND. &
5617	( lon_ful(ip) .LE. domaine_lon_max ) .AND. &
5618	( lat_ful(ip) .GE. domaine_lat_min ) .AND. &
5619	( lat_ful(ip) .LE. domaine_lat_max ) ) THEN
5620	!
5621	! look for point on GCM grid which this point on fine grid belongs to.
5622	! First look at the point on the model grid where we arrived just before. There is
5623	! a good chance that neighbouring points on the fine grid fall into the same model
5624	! grid box.
5625	!
5626	IF ( ( lon_ful(ip) .GE. lon_low(ilast) ) .AND. &
5627	( lon_ful(ip) .LT. lon_up(ilast) ) .AND. &
5628	( lat_ful(ip) .GE. lat_low(ilast) ) .AND. &
5629	( lat_ful(ip) .LT. lat_up(ilast) ) ) THEN
5630	!
5631	! We were lucky
5632	!
5633	IF (mask(ip) .GT. 0) THEN
5634	n_origpopden(ilast) = n_origpopden(ilast) + 1
5635	popd(ilast) = popd(ilast) + popden_orig(ip)
5636	ENDIF
5637	!
5638	ELSE
5639	!
5640	! Otherwise, look everywhere.
5641	! Begin close to last grid point.
5642	!
5643	found = .FALSE.
5644	idi = 1
5645	!
5646	DO WHILE ( (idi .LT. nbpt) .AND. ( .NOT. found ) )
5647
5648	!
5649	! forward and backward
5650	!
5651	DO ii = 1,2
5652	!
5653	IF ( ii .EQ. 1 ) THEN
5654	ib = ilast - idi
5655	ELSE
5656	ib = ilast + idi
5657	ENDIF
5658	!
5659	IF ( ( ib .GE. 1 ) .AND. ( ib .LE. nbpt ) ) THEN
5660	IF ( ( lon_ful(ip) .GE. lon_low(ib) ) .AND. &
5661	( lon_ful(ip) .LT. lon_up(ib) ) .AND. &
5662	( lat_ful(ip) .GE. lat_low(ib) ) .AND. &
5663	( lat_ful(ip) .LT. lat_up(ib) ) ) THEN
5664	!
5665	IF (mask(ip) .gt. 0) THEN
5666	popd(ib) = popd(ib) + popden_orig(ip)
5667	n_origpopden(ib) = n_origpopden(ib) + 1
5668	ENDIF
5669	ilast = ib
5670	found = .TRUE.
5671	!
5672	ENDIF
5673	ENDIF
5674	!
5675	ENDDO
5676	!
5677	idi = idi + 1
5678	!
5679	ENDDO
5680	!
5681	ENDIF ! lucky/not lucky
5682	!
5683	ENDIF ! in the domain
5684	ENDDO
5685
5686	WRITE(numout,*) ''
5687
5688	! determine fraction of popdening points in each box of the coarse grid
5689	DO ip=1,nbpt
5690	IF ( n_origpopden(ip) .GT. 0 ) THEN
5691	popd(ip) = popd(ip)/REAL(n_origpopden(ip),r_std)
5692	ELSE
5693	WRITE(numout,*) 'PROBLEM, no point in the popdens map found for this grid box'
5694	WRITE(numout,*) 'Longitude range : ', lon_low(ip), lon_up(ip)
5695	WRITE(numout,*) 'Latitude range : ', lat_low(ip), lat_up(ip)
5696	!
5697	WRITE(numout,*) 'Looking for nearest point on the popdens map file'
5698	CALL slowproc_nearest (ijml, lon_ful, lat_ful, &
5699	lalo(ip,2), lalo(ip,1), inear)
5700	WRITE(numout,*) 'Coordinates of the nearest point, ',inear,' :', &
5701	lon_ful(inear),lat_ful(inear)
5702	!
5703	popd(ip) = popden_orig(inear)
5704	ENDIF
5705	ENDDO
5706	!
5707	IF (printlev_loc>=1) WRITE(numout,*) '--',field_name,'--: Interpolation Done'
5708	WRITE(numout,*) ''
5709
5710	!
5711	DEALLOCATE(lon_up)
5712	DEALLOCATE(lon_low)
5713	DEALLOCATE(lat_up)
5714	DEALLOCATE(lat_low)
5715	DEALLOCATE(lat_ful)
5716	DEALLOCATE(lon_ful)
5717	DEALLOCATE(lat_lu)
5718	DEALLOCATE(lon_lu)
5719	DEALLOCATE(popden_lu)
5720	DEALLOCATE(popden_orig)
5721	DEALLOCATE(mask)
5722	!
5723	RETURN
5724	!
5725	END SUBROUTINE slowproc_read_annual
5726
5727	!endspit
5728
5729	END MODULE slowproc

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

Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_sechiba/slowproc.f90 @ 6737

Download in other formats: