Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_sechiba/slowproc.f90 @ 7599

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

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

Download in other formats: