Context Navigation

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

Last change on this file since 7709 was 7709, checked in by josefine.ghattas, 2 years ago

Integrated new irrigation scheme developed by Pedro Arboleda. See ticket #857
This corresponds to revsion 7708 of version pedro.arboleda/ORCHIDEE. Following differences were made but were not made on the pedro.arboleda/ORCHIDEE :

argumet place in call to routing_wrapper_intialize changed order
lines with only change in space were not taken
some indentation changed
set irrigation output as enalbled false if not do_irrigation

Property svn:keywords set to HeadURL Date Author Revision

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

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

Download in other formats: