Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90 @ 8418

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

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

Download in other formats: