Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

slowproc.f90

Last change on this file was 8579, checked in by antoine.bierjon, 2 weeks ago
Adapted read and interpolation of irrigation input files with XIOS_interpolation. Added a default value for AEI_SW interpolation (in case of missing value, with and without XIOS), which is chosen to be the average of the GMIA5 dataset from FAO. Related to ticket #857
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 248.9 KB

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

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

Context Navigation

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

Download in other formats: