Context Navigation

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

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

Reorganizing in output variables related to validation of interpolation, see ticket #812

Regrouped all output variables after interpolation into the new file sechiba_inter_diag
Renamed output variables axxx for availabilty for interpolation into interp_avail_axxx
Added some more variables needed to validate interpolation for pftmap and soils_param files
Set level to only output the final variables after interpolation
Removed variable "aveget5k" from field_def not sent from fortran code anymore

Property svn:keywords set to HeadURL Date Author Revision

File size: 215.2 KB

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

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

Download in other formats: