Context Navigation

source: branches/publications/ORCHIDEE_CAMEO_gmd_2022/src_sechiba/slowproc.f90 @ 8398

Last change on this file since 8398 was 7104, checked in by maureen.beaudor, 3 years ago
bug resolution : weight var in srf
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 267.1 KB

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

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

Download in other formats: