Context Navigation

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

Last change on this file since 8462 was 8462, checked in by josefine.ghattas, 4 months ago
The Moyano function describing the soil moisture effect on OM decomposition is added. It has been developed by Elodie Salmon in another branch and integrated in ORCHIDEE_2_2 by Bertrad Guenet. This commit corresponds to a corrected version of [8418].
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 245.4 KB

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

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

Download in other formats: