Context Navigation

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

Last change on this file since 7475 was 7475, checked in by josefine.ghattas, 2 years ago
Update orchidee.default using the tool create_orchidee_defalut.sh tool. Note that manual modifications in orchidee.default are overwritten. Reported some of the changes done previously in orchidee.default to the fortran source code. Note also that all sections in !Config must be there to have the tool working. Therfore added some missing !Config If lines.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 214.3 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_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	ELSE
1023	njsc = NINT(tmp_real)
1024	END IF
1025
1026	! Other soil parameters
1027
1028	var_name= 'ks'
1029	CALL ioconf_setatt_p('UNITS', 'mm/d')
1030	CALL ioconf_setatt_p('LONG_NAME','Soil saturated water content')
1031	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., ks, "gather", nbp_glo, index_g)
1032	IF ( ALL(ks(:) .EQ. val_exp ) ) THEN
1033	! ks is not in restart file
1034	call_slowproc_soilt=.TRUE.
1035	END IF
1036
1037	var_name= 'mcs'
1038	CALL ioconf_setatt_p('UNITS', 'm3/m3')
1039	CALL ioconf_setatt_p('LONG_NAME','')
1040	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcs, "gather", nbp_glo, index_g)
1041	IF ( ALL(mcs(:) .EQ. val_exp ) ) THEN
1042	! mcs is not in restart file
1043	call_slowproc_soilt=.TRUE.
1044	END IF
1045
1046	var_name= 'mcr'
1047	CALL ioconf_setatt_p('UNITS', 'm3/m3')
1048	CALL ioconf_setatt_p('LONG_NAME','')
1049	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcr, "gather", nbp_glo, index_g)
1050	IF ( ALL(mcr(:) .EQ. val_exp ) ) THEN
1051	! mcr is not in restart file
1052	call_slowproc_soilt=.TRUE.
1053	END IF
1054
1055	var_name= 'mcfc'
1056	CALL ioconf_setatt_p('UNITS', 'm3/m3')
1057	CALL ioconf_setatt_p('LONG_NAME','')
1058	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcfc, "gather", nbp_glo, index_g)
1059	IF ( ALL(mcfc(:) .EQ. val_exp ) ) THEN
1060	! mcfc is not in restart file
1061	call_slowproc_soilt=.TRUE.
1062	END IF
1063
1064	var_name= 'mcw'
1065	CALL ioconf_setatt_p('UNITS', 'm3/m3')
1066	CALL ioconf_setatt_p('LONG_NAME','')
1067	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcw, "gather", nbp_glo, index_g)
1068	IF ( ALL(mcw(:) .EQ. val_exp ) ) THEN
1069	! mcw is not in restart file
1070	call_slowproc_soilt=.TRUE.
1071	END IF
1072
1073	var_name= 'nvan'
1074	CALL ioconf_setatt_p('UNITS', '-')
1075	CALL ioconf_setatt_p('LONG_NAME','')
1076	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., nvan, "gather", nbp_glo, index_g)
1077	IF ( ALL(nvan(:) .EQ. val_exp ) ) THEN
1078	! nvan is not in restart file
1079	call_slowproc_soilt=.TRUE.
1080	END IF
1081
1082	var_name= 'avan'
1083	CALL ioconf_setatt_p('UNITS', 'm-1')
1084	CALL ioconf_setatt_p('LONG_NAME','')
1085	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., avan, "gather", nbp_glo, index_g)
1086	IF ( ALL(avan(:) .EQ. val_exp ) ) THEN
1087	! avan is not in restart file
1088	call_slowproc_soilt=.TRUE.
1089	END IF
1090
1091	var_name= 'clay_frac'
1092	CALL ioconf_setatt_p('UNITS', '-')
1093	CALL ioconf_setatt_p('LONG_NAME','Fraction of clay in each mesh')
1094	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., clayfraction, "gather", nbp_glo, index_g)
1095	IF ( ALL(clayfraction(:) .EQ. val_exp ) ) THEN
1096	! clayfraction is not in restart file
1097	call_slowproc_soilt=.TRUE.
1098	END IF
1099
1100	var_name= 'sand_frac'
1101	CALL ioconf_setatt_p('UNITS', '-')
1102	CALL ioconf_setatt_p('LONG_NAME','Fraction of sand in each mesh')
1103	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., sandfraction, "gather", nbp_glo, index_g)
1104	IF ( ALL(sandfraction(:) .EQ. val_exp ) ) THEN
1105	! sandfraction is not in restart file
1106	call_slowproc_soilt=.TRUE.
1107	END IF
1108
1109	! Calculate siltfraction not needed to be in restart file
1110	IF ( ALL( sandfraction(:) .EQ. val_exp) ) THEN
1111	siltfraction(:) = val_exp
1112	ELSE
1113	siltfraction(:) = 1. - clayfraction(:) - sandfraction(:)
1114	END IF
1115
1116	!! 2.b If we want to prescribe the soil parameters values with a unique value (e.g. site experiments)
1117	! AD: improve consistency with case spmipexp="unif"
1118
1119	IF (impsoilt) THEN
1120
1121	!Config Key = SOIL_FRACTIONS
1122	!Config Desc = Areal fraction of the 13 soil USDA textures; the dominant one is selected, Loam by default
1123	!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
1124	!Config If = IMPOSE_SOILT
1125	!Config Help = Determines the fraction of the 13 USDA classes with same order as in constantes_soil_var
1126	!Config Units = [-]
1127
1128	soilclass(:,:)=val_exp
1129	CALL setvar_p (soilclass, val_exp, 'SOIL_FRACTIONS', soilclass_default)
1130
1131	! Simplify a heterogeneous grid-cell into an homogeneous one
1132	! with the dominant texture
1133	njsc(:) = 0
1134	DO ji = 1, kjpindex
1135	! here we reduce to the dominant texture class
1136	njsc(ji) = MAXLOC(soilclass(ji,:),1)
1137	ENDDO
1138	njsc_imp = njsc(1) ! to prescribe the VG parameters consistely with the imposed texture
1139
1140	!Config Key = CLAY_FRACTION
1141	!Config Desc = Fraction of the clay fraction (0-dim mode)
1142	!Config Def = 0.2 if Loam
1143	!Config If = IMPOSE_SOIL
1144	!Config Help = Determines the fraction of clay in the grid box.
1145	!Config Units = [-]
1146	clayfraction_default = clayfrac_usda(njsc_imp)
1147	CALL setvar_p (clayfraction, val_exp, 'CLAY_FRACTION', clayfraction_default)
1148
1149	!Config Key = SAND_FRACTION
1150	!Config Desc = Fraction of the clay fraction (0-dim mode)
1151	!Config Def = 0.4
1152	!Config If = IMPOSE_SOIL
1153	!Config Help = Determines the fraction of clay in the grid box.
1154	!Config Units = [-]
1155	sandfraction_default = sandfrac_usda(njsc_imp)
1156	CALL setvar_p (sandfraction, val_exp, 'SAND_FRACTION', sandfraction_default)
1157
1158	! Calculate silt fraction
1159	siltfraction(:) = 1. - clayfraction(:) - sandfraction(:)
1160
1161	!Config Key = nvan
1162	!Config Desc = nvan parameter from Van genutchen equations
1163	!Config Def = 1.56 if Loam
1164	!Config If = IMPOSE_SOIL
1165	!Config Help = Determines the nvan in the grid box.
1166	!Config Units = [-]
1167	nvan_default = nvan_usda(njsc_imp)
1168	CALL setvar_p (nvan, val_exp, 'NVAN_IMP', nvan_default)
1169
1170	!Config Key = avan
1171	!Config Desc = avan parameter from Van genutchen equations
1172	!Config Def = 0.0036 if Loam
1173	!Config If = IMPOSE_SOIL
1174	!Config Help = Determines the avan in the grid box.
1175	!Config Units = [-]
1176	avan_default = avan_usda(njsc_imp)
1177	CALL setvar_p (avan, val_exp, 'AVAN_IMP', avan_default)
1178
1179	!Config Key = mcr
1180	!Config Desc = residual soil moisture
1181	!Config Def = 0.078 if Loam
1182	!Config If = IMPOSE_SOIL
1183	!Config Help = Determines the mcr in the grid box.
1184	!Config Units = [-]
1185	mcr_default = mcr_usda(njsc_imp)
1186	CALL setvar_p (mcr, val_exp, 'MCR_IMP', mcr_default)
1187
1188	!Config Key = mcs
1189	!Config Desc = saturation soil moisture
1190	!Config Def = 0.43 if Loam
1191	!Config If = IMPOSE_SOIL
1192	!Config Help = Determines the mcs in the grid box.
1193	!Config Units = [-]
1194	mcs_default = mcs_usda(njsc_imp)
1195	CALL setvar_p (mcs, val_exp, 'MCS_IMP', mcs_default)
1196
1197	!Config Key = ks
1198	!Config Desc = saturation conductivity
1199	!Config Def = 249.6 if Loam
1200	!Config If = IMPOSE_SOIL
1201	!Config Help = Determines the ks in the grid box.
1202	!Config Units = [mm/d]
1203	ks_default = ks_usda(njsc_imp)
1204	CALL setvar_p (ks, val_exp, 'KS_IMP', ks_default)
1205
1206	! By default, we calculate mcf and mcw from the above values, as in slowproc_soilt,
1207	! but they can be overruled by values from run.def
1208
1209	mvan(:) = un - (un / nvan(:))
1210	! Define matrix potential in mm for wilting point and field capacity (with sand vs clay-silt variation)
1211	psi_w(:) = 150000.
1212	DO ji=1, kjpindex
1213	IF ( ks(ji) .GE. 560 ) THEN ! Sandy soils (560 is equivalent of 2.75 at log scale of Ks, mm/d)
1214	psi_fc(ji) = 1000.
1215	ELSE ! Finer soils
1216	psi_fc(ji) = 3300.
1217	ENDIF
1218	ENDDO
1219	mcfc_default(:) = mcr(:) + (( mcs(:) - mcr(:)) / (un + ( avan(:) * psi_fc(:)) nvan(:)) mvan(:))
1220	mcw_default(:) = mcr(:) + (( mcs(:) - mcr(:)) / (un + ( avan(:) * psi_w(:)) nvan(:)) mvan(:))
1221
1222	!Config Key = mcfc
1223	!Config Desc = field capacity soil moisture
1224	!Config Def = 0.1654 if caclulated from default 5 parameters above
1225	!Config If = IMPOSE_SOIL
1226	!Config Help = Determines the mcfc in the grid box.
1227	!Config Units = [-]
1228	mcfc(:) = mcfc_default(:)
1229	CALL setvar_p (mcfc, val_exp, 'MCFC_IMP', mcfc_default)
1230
1231	!Config Key = mcw
1232	!Config Desc = wilting point soil moisture
1233	!Config Def = 0.0884 if caclulated from default 5 parameters above
1234	!Config If = IMPOSE_SOIL
1235	!Config Help = Determines the mcw in the grid box.
1236	!Config Units = [-]
1237	mcw(:) = mcw_default(:)
1238	CALL setvar_p (mcw, val_exp, 'MCW_IMP', mcw_default)
1239
1240	ELSE ! impsoilt = F
1241
1242	!! 2.c If some soil parameters are missing in the restart, we read a soil map
1243	! Note that slowproc_soilt calculates mcf and lcw as a function of avan, nvan, etc.
1244
1245	IF (call_slowproc_soilt) THEN
1246
1247	! If njsc is not in restart file, then initialize to default value
1248	! tbdone - reading from run.def file
1249	IF ( ALL(njsc(:) .EQ. undef_int )) THEN
1250	njsc(:) = usda_default ! 6 = Loam
1251	END IF
1252
1253	CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, &
1254	lalo, neighbours, resolution, contfrac, soilclass, &
1255	clayfraction, sandfraction, siltfraction)
1256
1257	call_slowproc_soilt=.FALSE.
1258
1259	ENDIF
1260
1261	ENDIF
1262
1263	! XIOS export of Ks before changing the vertical profile
1264	CALL xios_orchidee_send_field("ksref",ks) ! mm/d (for CMIP6, once)
1265
1266	!! 3. Read the infiltration related variables
1267	! This variable helps reducing surface runuff in flat areas
1268
1269	!! 3.a Looking first in the restart files
1270
1271	var_name= 'reinf_slope'
1272	CALL ioconf_setatt_p('UNITS', '-')
1273	CALL ioconf_setatt_p('LONG_NAME','Slope coef for reinfiltration')
1274	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinf_slope, "gather", nbp_glo, index_g)
1275
1276	!! 3.b We can also read from a map or prescribe, depending on GET_SLOPE
1277
1278	!Config Key = GET_SLOPE
1279	!Config Desc = Read slopes from file and do the interpolation
1280	!Config Def = n
1281	!Config If =
1282	!Config Help = Needed for reading the slopesfile and doing the interpolation. This will be
1283	! used by the re-infiltration parametrization
1284	!Config Units = [FLAG]
1285	get_slope = .FALSE.
1286	CALL getin_p('GET_SLOPE',get_slope)
1287
1288	IF (get_slope) THEN
1289
1290	!! Impose a constant value from the run.def (not anymore controlled by impveg)
1291
1292	!Config Key = REINF_SLOPE
1293	!Config Desc = Fraction of reinfiltrated surface runoff
1294	!Config Def = 0.1
1295	!Config If =
1296	!Config Help = Determines the reinfiltration ratio in the grid box due to flat areas
1297	!Config Units = [-]
1298	slope_default=0.1
1299	CALL setvar_p (reinf_slope, val_exp, 'SLOPE', slope_default)
1300
1301	ELSE
1302
1303	!! Else, if nof found from restart, we read from a map
1304	IF ( MINVAL(reinf_slope) .EQ. MAXVAL(reinf_slope) .AND. MAXVAL(reinf_slope) .EQ. val_exp .OR. get_slope) THEN
1305	IF (printlev_loc>=4) WRITE (numout,*) 'reinf_slope was not in restart file. Now call slowproc_slope'
1306
1307	CALL slowproc_slope(kjpindex, lalo, neighbours, resolution, contfrac, reinf_slope)
1308	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_slope'
1309
1310	ENDIF
1311
1312	ENDIF
1313
1314	!! 4. Read the vegetation related variables
1315
1316	!! 4.a General parameters
1317
1318	!Config Key = SECHIBA_QSINT
1319	!Config Desc = Interception reservoir coefficient
1320	!Config If = OK_SECHIBA
1321	!Config Def = 0.02
1322	!Config Help = Transforms leaf area index into size of interception reservoir
1323	!Config for slowproc_derivvar or stomate
1324	!Config Units = [m]
1325	CALL getin_p('SECHIBA_QSINT', qsintcst)
1326	IF (printlev >= 2) WRITE(numout, *)' SECHIBA_QSINT, qsintcst = ', qsintcst
1327
1328	!Config Key = SECHIBA_ZCANOP
1329	!Config Desc = Soil level used for canopy development (if STOMATE disactivated)
1330	!Config If = OK_SECHIBA and .NOT. OK_STOMATE
1331	!Config Def = 0.5
1332	!Config Help = The temperature at this soil depth is used to determine the LAI when
1333	!Config STOMATE is not activated.
1334	!Config Units = [m]
1335	zcanop = 0.5_r_std
1336	CALL setvar_p (zcanop, val_exp, 'SECHIBA_ZCANOP', 0.5_r_std)
1337
1338	! depth at center of the levels
1339	zsoil(1) = diaglev(1) / 2.
1340	DO l = 2, nslm
1341	zsoil(l) = ( diaglev(l) + diaglev(l-1) ) / 2.
1342	ENDDO
1343
1344	! index of this level
1345	vtmp = MINLOC ( ABS ( zcanop - zsoil(:) ) )
1346	lcanop = vtmp(1)
1347
1348	!! 4.b From restart first
1349
1350	found_restart=.TRUE.
1351	var_name= 'veget'
1352	CALL ioconf_setatt_p('UNITS', '-')
1353	CALL ioconf_setatt_p('LONG_NAME','Vegetation fraction')
1354	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget, "gather", nbp_glo, index_g)
1355	IF ( ALL( veget(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
1356
1357	var_name= 'veget_max'
1358	CALL ioconf_setatt_p('UNITS', '-')
1359	CALL ioconf_setatt_p('LONG_NAME','Maximum vegetation fraction')
1360	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., veget_max, "gather", nbp_glo, index_g)
1361	IF ( ALL( veget_max(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
1362
1363	IF (do_wood_harvest) THEN
1364	var_name= 'woodharvest'
1365	CALL ioconf_setatt_p('UNITS', 'gC m-2 yr-1')
1366	CALL ioconf_setatt_p('LONG_NAME','Harvest wood biomass')
1367	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., woodharvest, "gather", nbp_glo, index_g)
1368	IF ( ALL( woodharvest(:) .EQ. val_exp ) ) woodharvest(:)=zero
1369	END IF
1370
1371	var_name= 'frac_nobio'
1372	CALL ioconf_setatt_p('UNITS', '-')
1373	CALL ioconf_setatt_p('LONG_NAME','Special soil type fraction')
1374	CALL restget_p (rest_id, var_name, nbp_glo, nnobio, 1, kjit, .TRUE., frac_nobio, "gather", nbp_glo, index_g)
1375	IF ( ALL( frac_nobio(:,:) .EQ. val_exp ) ) found_restart=.FALSE.
1376
1377	IF (.NOT. impveg) THEN
1378	var_name= 'veget_year'
1379	CALL ioconf_setatt_p('UNITS', '-')
1380	CALL ioconf_setatt_p('LONG_NAME','Last year get in Land Use file.')
1381	IF (is_root_prc) THEN
1382	CALL restget (rest_id, var_name, 1 , 1 , 1, kjit, .TRUE., tmp_veget_year)
1383	IF (veget_reinit) THEN
1384	! Do not take the value read from restart file
1385	veget_year=veget_year_orig
1386	ELSE IF (tmp_veget_year(1) == val_exp) THEN
1387	! veget_year was not found in restart file
1388	veget_year=veget_year_orig
1389	ELSE
1390	! veget_year was found in restart file, transform to integer
1391	veget_year=INT(tmp_veget_year(1))
1392	ENDIF
1393	ENDIF
1394	CALL bcast(veget_year)
1395
1396	!
1397	!Config Key = VEGET_UPDATE
1398	!Config Desc = Update vegetation frequency: 0Y or 1Y
1399	!Config If =
1400	!Config Def = 0Y
1401	!Config Help = The veget datas will be update each this time step. Must be 0Y if IMPOSE_VEG=y.
1402	!Config Units = [years]
1403	!
1404	veget_update=0
1405	WRITE(veget_str,'(a)') '0Y'
1406	CALL getin_p('VEGET_UPDATE', veget_str)
1407	l=INDEX(TRIM(veget_str),'Y')
1408	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
1409	IF (printlev_loc >= 2) WRITE(numout,*) "Update frequency for land use in years :",veget_update
1410
1411	! Coherence test
1412	IF (veget_update > 0 .AND. ok_dgvm .AND. .NOT. agriculture) THEN
1413	CALL ipslerr_p(3,'slowproc_init',&
1414	'The combination DGVM=TRUE, AGRICULTURE=FALSE and VEGET_UPDATE>0 is not possible', &
1415	'Set VEGET_UPDATE=0Y in run.def','')
1416	END IF
1417
1418	ELSE
1419	! impveg=TRUE: there can not be any land use change, veget_update must be =0
1420	! Read VEGET_UPDATE from run.def and exit if it is different from 0Y
1421	veget_update=0
1422	WRITE(veget_str,'(a)') '0Y'
1423	CALL getin_p('VEGET_UPDATE', veget_str)
1424	l=INDEX(TRIM(veget_str),'Y')
1425	READ(veget_str(1:(l-1)),"(I2.2)") veget_update
1426	IF (veget_update /= 0) THEN
1427	WRITE(numout,*) 'veget_update=',veget_update,' is not coeherent with impveg=',impveg
1428	CALL ipslerr_p(3,'slowproc_init','Incoherent values between impveg and veget_update', &
1429	'VEGET_UPDATE must be equal to 0Y if IMPOSE_VEG=y (impveg=true)','')
1430	END IF
1431
1432	ENDIF
1433
1434	var_name= 'lai'
1435	CALL ioconf_setatt_p('UNITS', '-')
1436	CALL ioconf_setatt_p('LONG_NAME','Leaf area index')
1437	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., lai, "gather", nbp_glo, index_g)
1438
1439	! The height of the vegetation could in principle be recalculated at the beginning of the run.
1440	! However, this is very tedious, as many special cases have to be taken into account. This variable
1441	! is therefore saved in the restart file.
1442	var_name= 'height'
1443	CALL ioconf_setatt_p('UNITS', 'm')
1444	CALL ioconf_setatt_p('LONG_NAME','Height of vegetation')
1445	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., height, "gather", nbp_glo, index_g)
1446
1447	IF (read_lai)THEN
1448	var_name= 'laimap'
1449	CALL ioconf_setatt_p('UNITS', '-')
1450	CALL ioconf_setatt_p('LONG_NAME','Leaf area index read')
1451	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 12, kjit, .TRUE., laimap)
1452	ENDIF
1453
1454	CALL ioconf_setatt_p('UNITS', '-')
1455	CALL ioconf_setatt_p('LONG_NAME','Fraction of leaves in leaf age class ')
1456	DO jf = 1, nleafages
1457	! variable name is somewhat complicated as ioipsl does not allow 3d variables for the moment...
1458	WRITE(laistring,'(i4)') jf
1459	laistring=ADJUSTL(laistring)
1460	var_name='frac_age_'//laistring(1:LEN_TRIM(laistring))
1461	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE.,frac_age(:,:,jf), "gather", nbp_glo, index_g)
1462	ENDDO
1463
1464	!! 4.c Initialization of variables not found in restart file
1465
1466	IF ( impveg ) THEN
1467
1468	!! 4.1.a Case impveg=true: Initialization of variables by reading run.def
1469	!! The routine setvar_p will only initialize the variable if it was not found in restart file.
1470	!! We are on a point and thus we can read the information from the run.def
1471
1472	!Config Key = SECHIBA_VEGMAX
1473	!Config Desc = Maximum vegetation distribution within the mesh (0-dim mode)
1474	!Config If = IMPOSE_VEG
1475	!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
1476	!Config Help = The fraction of vegetation is read from the restart file. If
1477	!Config it is not found there we will use the values provided here.
1478	!Config Units = [-]
1479	CALL setvar_p (veget_max, val_exp, 'SECHIBA_VEGMAX', veget_ori_fixed_test_1)
1480
1481	!Config Key = SECHIBA_FRAC_NOBIO
1482	!Config Desc = Fraction of other surface types within the mesh (0-dim mode)
1483	!Config If = IMPOSE_VEG
1484	!Config Def = 0.0
1485	!Config Help = The fraction of ice, lakes, etc. is read from the restart file. If
1486	!Config it is not found there we will use the values provided here.
1487	!Config For the moment, there is only ice.
1488	!Config Units = [-]
1489	frac_nobio1 = frac_nobio(1,1)
1490	CALL setvar_p (frac_nobio1, val_exp, 'SECHIBA_FRAC_NOBIO', frac_nobio_fixed_test_1)
1491	frac_nobio(:,:) = frac_nobio1
1492
1493	IF (.NOT. found_restart) THEN
1494	! Call slowproc_veget to correct veget_max and to calculate veget and soiltiles
1495	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
1496	END IF
1497
1498	!Config Key = SECHIBA_LAI
1499	!Config Desc = LAI for all vegetation types (0-dim mode)
1500	!Config Def = 0., 8., 8., 4., 4.5, 4.5, 4., 4.5, 4., 2., 2., 2., 2.
1501	!Config If = IMPOSE_VEG
1502	!Config Help = The maximum LAI used in the 0dim mode. The values should be found
1503	!Config in the restart file. The new values of LAI will be computed anyway
1504	!Config at the end of the current day. The need for this variable is caused
1505	!Config by the fact that the model may stop during a day and thus we have not
1506	!Config yet been through the routines which compute the new surface conditions.
1507	!Config Units = [-]
1508	CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax)
1509
1510	!Config Key = SLOWPROC_HEIGHT
1511	!Config Desc = Height for all vegetation types
1512	!Config Def = 0., 30., 30., 20., 20., 20., 15., 15., 15., .5, .6, 1.0, 1.0
1513	!Config If = OK_SECHIBA
1514	!Config Help = The height used in the 0dim mode. The values should be found
1515	!Config in the restart file. The new values of height will be computed anyway
1516	!Config at the end of the current day. The need for this variable is caused
1517	!Config by the fact that the model may stop during a day and thus we have not
1518	!Config yet been through the routines which compute the new surface conditions.
1519	!Config Units = [m]
1520	CALL setvar_p (height, val_exp, 'SLOWPROC_HEIGHT', height_presc)
1521
1522
1523	ELSE IF ( .NOT. found_restart .OR. vegetmap_reset ) THEN
1524
1525	!! 4.1.b Case impveg=false and no restart files: Initialization by reading vegetation map
1526
1527	! Initialize veget_max and frac_nobio
1528	! Case without restart file
1529	IF (printlev_loc>=3) WRITE(numout,*) 'Before call slowproc_readvegetmax in initialization phase without restart files'
1530	IF (printlev_loc>=3) WRITE(numout,*) 'veget_year=', veget_year
1531
1532	! Call the routine to read the vegetation from file (output is veget_max_new)
1533	CALL slowproc_readvegetmax(kjpindex, lalo, neighbours, resolution, contfrac, &
1534	veget_max, veget_max_new, frac_nobio_new, veget_year, .TRUE.)
1535	IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_readvegetmax in initialization phase'
1536
1537	! Update vegetation with values read from the file
1538	veget_max = veget_max_new
1539	frac_nobio = frac_nobio_new
1540
1541	IF (do_wood_harvest) THEN
1542	! Read the new the wood harvest map from file. Output is wood harvest
1543	CALL slowproc_woodharvest(kjpindex, lalo, neighbours, resolution, contfrac, woodharvest)
1544	ENDIF
1545
1546	!! Reset totaly or partialy veget_max if using DGVM
1547	IF ( ok_dgvm ) THEN
1548	! If we are dealing with dynamic vegetation then all natural PFTs should be set to veget_max = 0
1549	! In case no agriculture is desired, agriculture PFTS should be set to 0 as well
1550	IF (agriculture) THEN
1551	DO jv = 2, nvm
1552	IF (natural(jv)) THEN
1553	veget_max(:,jv)=zero
1554	ENDIF
1555	ENDDO
1556
1557	! Calculate the fraction of crop for each point.
1558	! Sum only on the indexes corresponding to the non_natural pfts
1559	frac_crop_tot(:) = zero
1560	DO jv = 2, nvm
1561	IF(.NOT. natural(jv)) THEN
1562	DO ji = 1, kjpindex
1563	frac_crop_tot(ji) = frac_crop_tot(ji) + veget_max(ji,jv)
1564	ENDDO
1565	ENDIF
1566	END DO
1567
1568	! Calculate the fraction of bare soil
1569	DO ji = 1, kjpindex
1570	veget_max(ji,1) = un - frac_crop_tot(ji) - SUM(frac_nobio(ji,:))
1571	ENDDO
1572	ELSE
1573	veget_max(:,:) = zero
1574	DO ji = 1, kjpindex
1575	veget_max(ji,1) = un - SUM(frac_nobio(ji,:))
1576	ENDDO
1577	END IF
1578	END IF ! end ok_dgvm
1579
1580
1581	! Call slowproc_veget to correct veget_max and to calculate veget and soiltiles
1582	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
1583
1584	END IF ! end impveg
1585
1586	!! 4.d Continue initializing variables not found in restart file. Case for both impveg=true and false.
1587
1588	! Initialize laimap for the case read_lai if not found in restart file
1589	IF (read_lai) THEN
1590	IF ( ALL( laimap(:,:,:) .EQ. val_exp) ) THEN
1591	! Interpolation of LAI
1592	CALL slowproc_interlai (kjpindex, lalo, resolution, neighbours, contfrac, laimap)
1593	ENDIF
1594	ENDIF
1595
1596	! Initialize lai if not found in restart file and not already initialized using impveg
1597	IF ( MINVAL(lai) .EQ. MAXVAL(lai) .AND. MAXVAL(lai) .EQ. val_exp) THEN
1598	IF (read_lai) THEN
1599	stempdiag2_bid(1:kjpindex,1:nslm) = stempdiag_bid
1600	CALL slowproc_lai (kjpindex, lcanop, stempdiag2_bid, &
1601	lalo,resolution,lai,laimap)
1602	ELSE
1603	! If we start from scratch, we set lai to zero for consistency with stomate
1604	lai(:,:) = zero
1605	ENDIF
1606
1607	frac_age(:,:,1) = un
1608	frac_age(:,:,2) = zero
1609	frac_age(:,:,3) = zero
1610	frac_age(:,:,4) = zero
1611	ENDIF
1612
1613	! Initialize heigth if not found in restart file and not already initialized using impveg
1614	IF ( MINVAL(height) .EQ. MAXVAL(height) .AND. MAXVAL(height) .EQ. val_exp) THEN
1615	! Impose height
1616	DO jv = 1, nvm
1617	height(:,jv) = height_presc(jv)
1618	ENDDO
1619	ENDIF
1620
1621	!! 5. Some calculations always done, with and without restart files
1622
1623	! The variables veget, veget_max and frac_nobio were all read from restart file or initialized above.
1624	! Calculate now totfrac_nobio and soiltiles using these variables.
1625
1626	! Calculate totfrac_nobio
1627	totfrac_nobio(:) = zero
1628	DO jv = 1, nnobio
1629	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
1630	ENDDO
1631
1632	! Calculate soiltile. This variable do not need to be in the restart file.
1633	! The sum of all soiltiles makes one, and corresponds to the bio fraction
1634	! of the grid cell (called vegtot in hydrol)
1635	soiltile(:,:) = zero
1636	DO jv = 1, nvm
1637	jst = pref_soil_veg(jv)
1638	DO ji = 1, kjpindex
1639	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
1640	ENDDO
1641	ENDDO
1642	DO ji = 1, kjpindex
1643	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
1644	soiltile(ji,:)=soiltile(ji,:)/(1-totfrac_nobio(ji))
1645	ENDIF
1646	ENDDO
1647
1648	! Always calculate tot_bare_soil
1649	! Fraction of bare soil in the mesh (bio+nobio)
1650	tot_bare_soil(:) = veget_max(:,1)
1651	DO jv = 2, nvm
1652	DO ji =1, kjpindex
1653	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
1654	ENDDO
1655	END DO
1656
1657
1658	!! Calculate fraction of landuse tiles to be used only for diagnostic variables
1659	fraclut(:,:)=0
1660	nwdFraclut(:,id_psl)=0
1661	nwdFraclut(:,id_crp)=1.
1662	nwdFraclut(:,id_urb)=xios_default_val
1663	nwdFraclut(:,id_pst)=xios_default_val
1664	DO jv=1,nvm
1665	IF (natural(jv)) THEN
1666	fraclut(:,id_psl) = fraclut(:,id_psl) + veget_max(:,jv)
1667	IF(.NOT. is_tree(jv)) THEN
1668	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl) + veget_max(:,jv)
1669	ENDIF
1670	ELSE
1671	fraclut(:,id_crp) = fraclut(:,id_crp) + veget_max(:,jv)
1672	ENDIF
1673	END DO
1674
1675	WHERE (fraclut(:,id_psl) > min_sechiba)
1676	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl)/fraclut(:,id_psl)
1677	ELSEWHERE
1678	nwdFraclut(:,id_psl) = xios_default_val
1679	END WHERE
1680
1681
1682	IF (printlev_loc>=3) WRITE (numout,*) ' slowproc_init done '
1683
1684	END SUBROUTINE slowproc_init
1685
1686	!! ================================================================================================================================
1687	!! SUBROUTINE : slowproc_clear
1688	!!
1689	!>\BRIEF Clear all variables related to slowproc and stomate modules
1690	!!
1691	!_ ================================================================================================================================
1692
1693	SUBROUTINE slowproc_clear
1694
1695	! 1 clear all the variables defined as common for the routines in slowproc
1696
1697	IF (ALLOCATED (clayfraction)) DEALLOCATE (clayfraction)
1698	IF (ALLOCATED (sandfraction)) DEALLOCATE (sandfraction)
1699	IF (ALLOCATED (siltfraction)) DEALLOCATE (siltfraction)
1700	IF (ALLOCATED (laimap)) DEALLOCATE (laimap)
1701	IF (ALLOCATED (veget_max_new)) DEALLOCATE (veget_max_new)
1702	IF (ALLOCATED (woodharvest)) DEALLOCATE (woodharvest)
1703	IF (ALLOCATED (frac_nobio_new)) DEALLOCATE (frac_nobio_new)
1704
1705	! 2. Clear all the variables in stomate
1706
1707	CALL stomate_clear
1708	!
1709	END SUBROUTINE slowproc_clear
1710
1711	!! ================================================================================================================================
1712	!! SUBROUTINE : slowproc_derivvar
1713	!!
1714	!>\BRIEF Initializes variables related to the
1715	!! parameters to be assimilated, the maximum water on vegetation, the vegetation height,
1716	!! and the fraction of soil covered by dead leaves and the vegetation height
1717	!!
1718	!! DESCRIPTION : (definitions, functional, design, flags):
1719	!! (1) Initialization of the variables relevant for the assimilation parameters
1720	!! (2) Intialization of the fraction of soil covered by dead leaves
1721	!! (3) Initialization of the Vegetation height per PFT
1722	!! (3) Initialization the maximum water on vegetation for interception with a particular treatement of the PFT no.1
1723	!!
1724	!! RECENT CHANGE(S): None
1725	!!
1726	!! MAIN OUTPUT VARIABLE(S): ::qsintmax, ::deadleaf_cover, ::assim_param, ::height
1727	!!
1728	!! REFERENCE(S) : None
1729	!!
1730	!! FLOWCHART : None
1731	!! \n
1732	!_ ================================================================================================================================
1733
1734	SUBROUTINE slowproc_derivvar (kjpindex, veget, lai, &
1735	qsintmax, deadleaf_cover, assim_param, height, temp_growth)
1736
1737	!! INTERFACE DESCRIPTION
1738
1739	!! 0.1 Input scalar and fields
1740	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size - terrestrial pixels only
1741	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
1742	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: lai !! PFT leaf area index (m^{2} m^{-2})
1743
1744	!! 0.2. Output scalar and fields
1745	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintmax !! Maximum water on vegetation for interception(mm)
1746	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: deadleaf_cover !! fraction of soil covered by dead leaves (unitless)
1747	REAL(r_std), DIMENSION (kjpindex,nvm,npco2), INTENT (out) :: assim_param !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
1748	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: height !! height of the vegetation or surface in general ??? (m)
1749	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_growth !! growth temperature (Â°C)
1750	!
1751	!! 0.3 Local declaration
1752	INTEGER(i_std) :: jv !! Local indices
1753	!_ ================================================================================================================================
1754
1755	!
1756	! 1. Initialize (why here ??) the variables revelant for the assimilation parameters
1757	!
1758	DO jv = 1, nvm
1759	assim_param(:,jv,ivcmax) = vcmax_fix(jv)
1760	ENDDO
1761
1762	!
1763	! 2. Intialize the fraction of soil covered by dead leaves
1764	!
1765	deadleaf_cover(:) = zero
1766
1767	!
1768	! 3. Initialize the Vegetation height per PFT
1769	!
1770	DO jv = 1, nvm
1771	height(:,jv) = height_presc(jv)
1772	ENDDO
1773	!
1774	! 4. Initialize the maximum water on vegetation for interception
1775	!
1776	qsintmax(:,:) = qsintcst * veget(:,:) * lai(:,:)
1777
1778	! Added by Nathalie - July 2006
1779	! Initialize the case of the PFT no.1 to zero
1780	qsintmax(:,1) = zero
1781
1782	temp_growth(:)=25.
1783
1784	END SUBROUTINE slowproc_derivvar
1785
1786
1787	!! ================================================================================================================================
1788	!! SUBROUTINE : slowproc_mean
1789	!!
1790	!>\BRIEF Accumulates field_in over a period of dt_tot.
1791	!! Has to be called at every time step (dt).
1792	!! Mean value is calculated if ldmean=.TRUE.
1793	!! field_mean must be initialized outside of this routine!
1794	!!
1795	!! DESCRIPTION : (definitions, functional, design, flags):
1796	!! (1) AcumAcuumlm
1797	!!
1798	!! RECENT CHANGE(S): None
1799	!!
1800	!! MAIN OUTPUT VARIABLE(S): ::field_main
1801	!!
1802	!! REFERENCE(S) : None
1803	!!
1804	!! FLOWCHART : None
1805	!! \n
1806	!_ ================================================================================================================================
1807
1808	SUBROUTINE slowproc_mean (npts, n_dim2, dt_tot, dt, ldmean, field_in, field_mean)
1809
1810	!
1811	!! 0 declarations
1812
1813	!! 0.1 input scalar and variables
1814	INTEGER(i_std), INTENT(in) :: npts !! Domain size- terrestrial pixels only
1815	INTEGER(i_std), INTENT(in) :: n_dim2 !! Number of PFTs
1816	REAL(r_std), INTENT(in) :: dt_tot !! Time step of stomate (in days). The period over which the accumulation or the mean is computed
1817	REAL(r_std), INTENT(in) :: dt !! Time step in days
1818	LOGICAL, INTENT(in) :: ldmean !! Flag to calculate the mean after the accumulation ???
1819	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_in !! Daily field
1820
1821	!! 0.3 Modified field; The computed sum or mean field over dt_tot time period depending on the flag ldmean
1822	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_mean !! Accumulated field at dt_tot time period or mean field over dt_tot
1823
1824
1825	!_ ================================================================================================================================
1826
1827	!
1828	! 1. Accumulation the field over dt_tot period
1829	!
1830	field_mean(:,:) = field_mean(:,:) + field_in(:,:) * dt
1831
1832	!
1833	! 2. If the flag ldmean set, the mean field is computed over dt_tot period
1834	!
1835	IF (ldmean) THEN
1836	field_mean(:,:) = field_mean(:,:) / dt_tot
1837	ENDIF
1838
1839	END SUBROUTINE slowproc_mean
1840
1841
1842
1843	!! ================================================================================================================================
1844	!! SUBROUTINE : slowproc_long
1845	!!
1846	!>\BRIEF Calculates a temporally smoothed field (field_long) from
1847	!! instantaneous input fields.Time constant tau determines the strength of the smoothing.
1848	!! For tau -> infinity??, field_long becomes the true mean value of field_inst
1849	!! (but the spinup becomes infinietly long, too).
1850	!! field_long must be initialized outside of this routine!
1851	!!
1852	!! DESCRIPTION : (definitions, functional, design, flags):
1853	!! (1) Testing the time coherence betwen the time step dt and the time tau over which
1854	!! the rescaled of the mean is performed
1855	!! (2) Computing the rescaled mean over tau period
1856	!! MAIN OUTPUT VARIABLE(S): field_long
1857	!!
1858	!! RECENT CHANGE(S): None
1859	!!
1860	!! MAIN OUTPUT VARIABLE(S): ::field_long
1861	!!
1862	!! REFERENCE(S) : None
1863	!!
1864	!! FLOWCHART : None
1865	!! \n
1866	!_ ================================================================================================================================
1867
1868	SUBROUTINE slowproc_long (npts, n_dim2, dt, tau, field_inst, field_long)
1869
1870	!
1871	! 0 declarations
1872	!
1873
1874	! 0.1 input scalar and fields
1875
1876	INTEGER(i_std), INTENT(in) :: npts !! Domain size- terrestrial pixels only
1877	INTEGER(i_std), INTENT(in) :: n_dim2 !! Second dimension of the fields, which represents the number of PFTs
1878	REAL(r_std), INTENT(in) :: dt !! Time step in days
1879	REAL(r_std), INTENT(in) :: tau !! Integration time constant (has to have same unit as dt!)
1880	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(in) :: field_inst !! Instantaneous field
1881
1882
1883	! 0.2 modified field
1884
1885	! Long-term field
1886	REAL(r_std), DIMENSION(npts,n_dim2), INTENT(inout) :: field_long !! Mean value of the instantaneous field rescaled at tau time period
1887
1888	!_ ================================================================================================================================
1889
1890	!
1891	! 1 test coherence of the time
1892
1893	IF ( ( tau .LT. dt ) .OR. ( dt .LE. zero ) .OR. ( tau .LE. zero ) ) THEN
1894	WRITE(numout,*) 'slowproc_long: Problem with time steps'
1895	WRITE(numout,*) 'dt=',dt
1896	WRITE(numout,*) 'tau=',tau
1897	ENDIF
1898
1899	!
1900	! 2 integration of the field over tau
1901
1902	field_long(:,:) = ( field_inst(:,:)dt + field_long(:,:)(tau-dt) ) / tau
1903
1904	END SUBROUTINE slowproc_long
1905
1906
1907	!! ================================================================================================================================
1908	!! SUBROUTINE : slowproc_veget_max_limit
1909	!!
1910	!>\BRIEF Set small fractions of veget_max to zero and normalize to keep the sum equal 1
1911	!!
1912	!! DESCRIPTION : Set small fractions of veget_max to zero and normalize to keep the sum equal 1
1913	!!
1914	!! RECENT CHANGE(S): The subroutine was previously a part of slowproc_veget,
1915	!! but was separated to be called also from slowproc_readvegetmax in order
1916	!! to have limited/normalized vegetation fractions right after its reading
1917	!! from the file (added by V.Bastrikov, 15/06/2019)
1918	!!
1919	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, veget_max
1920	!!
1921	!! REFERENCE(S) : None
1922	!!
1923	!! FLOWCHART : None
1924	!! \n
1925	!_ ================================================================================================================================
1926
1927	SUBROUTINE slowproc_veget_max_limit (kjpindex, frac_nobio, veget_max)
1928	!
1929	! 0. Declarations
1930	!
1931	! 0.1 Input variables
1932	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1933
1934	! 0.2 Modified variables
1935	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
1936	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
1937
1938	! 0.4 Local scalar and varaiables
1939	INTEGER(i_std) :: ji, jv !! indices
1940	REAL(r_std) :: SUMveg !! Total vegetation summed across PFTs
1941
1942	!_ ================================================================================================================================
1943	IF (printlev_loc >= 3) WRITE(numout,*) 'Entering slowproc_veget_max_limit'
1944
1945	!! Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
1946	DO ji = 1, kjpindex
1947	IF ( SUM(frac_nobio(ji,:)) .LT. min_vegfrac ) THEN
1948	frac_nobio(ji,:) = zero
1949	ENDIF
1950
1951	IF (.NOT. ok_dgvm) THEN
1952	DO jv = 1, nvm
1953	IF ( veget_max(ji,jv) .LT. min_vegfrac ) THEN
1954	veget_max(ji,jv) = zero
1955	ENDIF
1956	ENDDO
1957	END IF
1958
1959	!! Normalize to keep the sum equal 1.
1960	SUMveg = SUM(frac_nobio(ji,:))+SUM(veget_max(ji,:))
1961	frac_nobio(ji,:) = frac_nobio(ji,:)/SUMveg
1962	veget_max(ji,:) = veget_max(ji,:)/SUMveg
1963	ENDDO
1964
1965	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_veget_max_limit ended'
1966
1967	END SUBROUTINE slowproc_veget_max_limit
1968
1969
1970	!! ================================================================================================================================
1971	!! SUBROUTINE : slowproc_veget
1972	!!
1973	!>\BRIEF Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
1974	!!
1975	!! DESCRIPTION : Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
1976	!! (1) Set veget_max and frac_nobio for fraction smaller than min_vegfrac.
1977	!! (2) Calculate veget
1978	!! (3) Calculate totfrac_nobio
1979	!! (4) Calculate soiltile
1980	!! (5) Calculate fraclut
1981	!!
1982	!! RECENT CHANGE(S): None
1983	!!
1984	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut
1985	!!
1986	!! REFERENCE(S) : None
1987	!!
1988	!! FLOWCHART : None
1989	!! \n
1990	!_ ================================================================================================================================
1991
1992	SUBROUTINE slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
1993	!
1994	! 0. Declarations
1995	!
1996	! 0.1 Input variables
1997	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1998	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! PFT leaf area index (m^{2} m^{-2})
1999
2000	! 0.2 Modified variables
2001	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
2002	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
2003
2004	! 0.3 Output variables
2005	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
2006	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio
2007	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
2008	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
2009	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
2010
2011	! 0.4 Local scalar and varaiables
2012	INTEGER(i_std) :: ji, jv, jst !! indices
2013
2014	!_ ================================================================================================================================
2015	IF (printlev_loc > 8) WRITE(numout,*) 'Entering slowproc_veget'
2016
2017	!! 1. Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
2018	!! Normalize to have the sum equal 1.
2019	CALL slowproc_veget_max_limit(kjpindex, frac_nobio, veget_max)
2020
2021	!! 2. Calculate veget
2022	!! If lai of a vegetation type (jv > 1) is small, increase soil part
2023	!! stomate-like calculation
2024	DO ji = 1, kjpindex
2025	veget(ji,1)=veget_max(ji,1)
2026	DO jv = 2, nvm
2027	veget(ji,jv) = veget_max(ji,jv) * ( un - exp( - lai(ji,jv) * ext_coeff_vegetfrac(jv) ) )
2028	ENDDO
2029	ENDDO
2030
2031
2032	!! 3. Calculate totfrac_nobio
2033	totfrac_nobio(:) = zero
2034	DO jv = 1, nnobio
2035	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
2036	ENDDO
2037
2038
2039	!! 4. Calculate soiltiles
2040	!! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
2041	!! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
2042	!! in the other modules to perform separated energy balances
2043	! The sum of all soiltiles makes one, and corresponds to the bio fraction
2044	! of the grid cell (called vegtot in hydrol)
2045	soiltile(:,:) = zero
2046	DO jv = 1, nvm
2047	jst = pref_soil_veg(jv)
2048	DO ji = 1, kjpindex
2049	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
2050	ENDDO
2051	ENDDO
2052	DO ji = 1, kjpindex
2053	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
2054	soiltile(ji,:)=soiltile(ji,:)/(1.-totfrac_nobio(ji))
2055	ENDIF
2056	ENDDO
2057
2058	!! 5. Calculate fraction of landuse tiles to be used only for diagnostic variables
2059	fraclut(:,:)=0
2060	nwdFraclut(:,id_psl)=0
2061	nwdFraclut(:,id_crp)=1.
2062	nwdFraclut(:,id_urb)=xios_default_val
2063	nwdFraclut(:,id_pst)=xios_default_val
2064	DO jv=1,nvm
2065	IF (natural(jv)) THEN
2066	fraclut(:,id_psl) = fraclut(:,id_psl) + veget_max(:,jv)
2067	IF(.NOT. is_tree(jv)) THEN
2068	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl) + veget_max(:,jv)
2069	ENDIF
2070	ELSE
2071	fraclut(:,id_crp) = fraclut(:,id_crp) + veget_max(:,jv)
2072	ENDIF
2073	END DO
2074
2075	WHERE (fraclut(:,id_psl) > min_sechiba)
2076	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl)/fraclut(:,id_psl)
2077	ELSEWHERE
2078	nwdFraclut(:,id_psl) = xios_default_val
2079	END WHERE
2080
2081	END SUBROUTINE slowproc_veget
2082
2083
2084	!! ================================================================================================================================
2085	!! SUBROUTINE : slowproc_lai
2086	!!
2087	!>\BRIEF Do the interpolation of lai for the PFTs in case the laimap is not read
2088	!!
2089	!! DESCRIPTION : (definitions, functional, design, flags):
2090	!! (1) Interplation by using the mean value of laimin and laimax for the PFTs
2091	!! (2) Interpolation between laimax and laimin values by using the temporal
2092	!! variations
2093	!! (3) If problem occurs during the interpolation, the routine stops
2094	!!
2095	!! RECENT CHANGE(S): None
2096	!!
2097	!! MAIN OUTPUT VARIABLE(S): ::lai
2098	!!
2099	!! REFERENCE(S) : None
2100	!!
2101	!! FLOWCHART : None
2102	!! \n
2103	!_ ================================================================================================================================
2104
2105	SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,laimap)
2106	!
2107	! 0. Declarations
2108	!
2109	!! 0.1 Input variables
2110	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
2111	INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI
2112	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag !! Soil temperature (K) ???
2113	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
2114	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x an y of the grid (m) - surface area of the gridbox
2115	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: laimap !! map of lai read
2116
2117	!! 0.2 Output
2118	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! PFT leaf area index (m^{2} m^{-2})LAI
2119
2120	!! 0.4 Local
2121	INTEGER(i_std) :: ji,jv !! Local indices
2122	!_ ================================================================================================================================
2123
2124	!
2125	IF ( .NOT. read_lai ) THEN
2126
2127	lai(: ,1) = zero
2128	! On boucle sur 2,nvm au lieu de 1,nvm
2129	DO jv = 2,nvm
2130	SELECT CASE (type_of_lai(jv))
2131
2132	CASE ("mean ")
2133	!
2134	! 1. do the interpolation between laimax and laimin
2135	!
2136	lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv))
2137	!
2138	CASE ("inter")
2139	!
2140	! 2. do the interpolation between laimax and laimin
2141	!
2142	DO ji = 1,kjpindex
2143	lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv))
2144	ENDDO
2145	!
2146	CASE default
2147	!
2148	! 3. Problem
2149	!
2150	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
2151	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
2152	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=false','')
2153	END SELECT
2154
2155	ENDDO
2156	!
2157	ELSE
2158	lai(: ,1) = zero
2159	! On boucle sur 2,nvm au lieu de 1,nvm
2160	DO jv = 2,nvm
2161
2162	SELECT CASE (type_of_lai(jv))
2163
2164	CASE ("mean ")
2165	!
2166	! 1. force MAXVAL of laimap on lai on this PFT
2167	!
2168	DO ji = 1,kjpindex
2169	lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
2170	ENDDO
2171	!
2172	CASE ("inter")
2173	!
2174	! 2. do the interpolation between laimax and laimin
2175	!
2176	!
2177	! If January
2178	!
2179	IF (month_end .EQ. 1 ) THEN
2180	IF (day_end .LE. 15) THEN
2181	lai(:,jv) = laimap(:,jv,12)(1-(day_end+15)/30.) + laimap(:,jv,1)((day_end+15)/30.)
2182	ELSE
2183	lai(:,jv) = laimap(:,jv,1)(1-(day_end-15)/30.) + laimap(:,jv,2)((day_end-15)/30.)
2184	ENDIF
2185	!
2186	! If December
2187	!
2188	ELSE IF (month_end .EQ. 12) THEN
2189	IF (day_end .LE. 15) THEN
2190	lai(:,jv) = laimap(:,jv,11)(1-(day_end+15)/30.) + laimap(:,jv,12)((day_end+15)/30.)
2191	ELSE
2192	lai(:,jv) = laimap(:,jv,12)(1-(day_end-15)/30.) + laimap(:,jv,1)((day_end-15)/30.)
2193	ENDIF
2194	!
2195	! ELSE
2196	!
2197	ELSE
2198	IF (day_end .LE. 15) THEN
2199	lai(:,jv) = laimap(:,jv,month_end-1)(1-(day_end+15)/30.) + laimap(:,jv,month_end)((day_end+15)/30.)
2200	ELSE
2201	lai(:,jv) = laimap(:,jv,month_end)(1-(day_end-15)/30.) + laimap(:,jv,month_end+1)((day_end-15)/30.)
2202	ENDIF
2203	ENDIF
2204	!
2205	CASE default
2206	!
2207	! 3. Problem
2208	!
2209	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
2210	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
2211	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=true','')
2212	END SELECT
2213
2214	ENDDO
2215	ENDIF
2216
2217	END SUBROUTINE slowproc_lai
2218
2219	!! ================================================================================================================================
2220	!! SUBROUTINE : slowproc_interlai
2221	!!
2222	!>\BRIEF Interpolate the LAI map to the grid of the model
2223	!!
2224	!! DESCRIPTION : (definitions, functional, design, flags):
2225	!!
2226	!! RECENT CHANGE(S): None
2227	!!
2228	!! MAIN OUTPUT VARIABLE(S): ::laimap
2229	!!
2230	!! REFERENCE(S) : None
2231	!!
2232	!! FLOWCHART : None
2233	!! \n
2234	!_ ================================================================================================================================
2235
2236	SUBROUTINE slowproc_interlai(nbpt, lalo, resolution, neighbours, contfrac, laimap)
2237
2238	USE interpweight
2239
2240	IMPLICIT NONE
2241
2242	!
2243	!
2244	!
2245	! 0.1 INPUT
2246	!
2247	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
2248	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
2249	!! (beware of the order = 1 : latitude, 2 : longitude)
2250	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
2251	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2252	!! (1=North and then clockwise)
2253	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
2254	!
2255	! 0.2 OUTPUT
2256	!
2257	REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) !! lai read variable and re-dimensioned
2258	!
2259	! 0.3 LOCAL
2260	!
2261	CHARACTER(LEN=80) :: filename !! name of the LAI map read
2262	INTEGER(i_std) :: ib, ip, jp, it, jv
2263	REAL(r_std) :: lmax, lmin, ldelta
2264	LOGICAL :: renormelize_lai ! flag to force LAI renormelization
2265	INTEGER :: ier
2266
2267	REAL(r_std), DIMENSION(nbpt) :: alaimap !! availability of the lai interpolation
2268	INTEGER, DIMENSION(4) :: invardims
2269	REAL(r_std), DIMENSION(nbpt,nvm,12) :: lairefrac !! lai fractions re-dimensioned
2270	REAL(r_std), DIMENSION(nbpt,nvm,12) :: fraclaiinterp !! lai fractions re-dimensioned
2271	REAL(r_std), DIMENSION(:), ALLOCATABLE :: vmin, vmax !! min/max values to use for the
2272	!! renormalization
2273	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2274	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2275	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
2276	!! (variabletypevals(1) = -un, not used)
2277	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2278	!! 'XYKindTime': Input values are kinds
2279	!! of something with a temporal
2280	!! evolution on the dx*dy matrix'
2281	LOGICAL :: nonegative !! whether negative values should be removed
2282	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2283	!! 'nomask': no-mask is applied
2284	!! 'mbelow': take values below maskvals(1)
2285	!! 'mabove': take values above maskvals(1)
2286	!! 'msumrange': take values within 2 ranges;
2287	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2288	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2289	!! (normalized by maskvals(3))
2290	!! 'var': mask values are taken from a
2291	!! variable inside the file (>0)
2292	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2293	!! `maskingtype')
2294	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2295	!_ ================================================================================================================================
2296
2297	!
2298	!Config Key = LAI_FILE
2299	!Config Desc = Name of file from which the vegetation map is to be read
2300	!Config If = LAI_MAP
2301	!Config Def = lai2D.nc
2302	!Config Help = The name of the file to be opened to read the LAI
2303	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
2304	!Config map which is derived from a Nicolas VIOVY one.
2305	!Config Units = [FILE]
2306	!
2307	filename = 'lai2D.nc'
2308	CALL getin_p('LAI_FILE',filename)
2309	variablename = 'LAI'
2310
2311	IF (xios_interpolation) THEN
2312	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_interlai: Use XIOS to read and interpolate " &
2313	// TRIM(filename) //" for variable " //TRIM(variablename)
2314
2315	CALL xios_orchidee_recv_field('lai_interp',lairefrac)
2316	CALL xios_orchidee_recv_field('frac_lai_interp',fraclaiinterp)
2317	alaimap(:) = fraclaiinterp(:,1,1)
2318	ELSE
2319
2320	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_interlai: Start interpolate " &
2321	// TRIM(filename) //" for variable " //TRIM(variablename)
2322
2323	! invardims: shape of variable in input file to interpolate
2324	invardims = interpweight_get_var4dims_file(filename, variablename)
2325	! Check coherence of dimensions read from the file
2326	IF (invardims(4) /= 12) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of time dimension in input file for lai','','')
2327	IF (invardims(3) /= nvm) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of PFT dimension in input file for lai','','')
2328
2329	ALLOCATE(vmin(nvm),stat=ier)
2330	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmin','','')
2331
2332	ALLOCATE(vmax(nvm), STAT=ier)
2333	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmax','','')
2334
2335
2336	! Assigning values to vmin, vmax
2337	vmin = un
2338	vmax = nvm*un
2339
2340	variabletypevals = -un
2341
2342	!! Variables for interpweight
2343	! Type of calculation of cell fractions
2344	fractype = 'default'
2345	! Name of the longitude and latitude in the input file
2346	lonname = 'longitude'
2347	latname = 'latitude'
2348	! Should negative values be set to zero from input file?
2349	nonegative = .TRUE.
2350	! Type of mask to apply to the input data (see header for more details)
2351	maskingtype = 'mbelow'
2352	! Values to use for the masking
2353	maskvals = (/ 20., undef_sechiba, undef_sechiba /)
2354	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2355	namemaskvar = ''
2356
2357	CALL interpweight_4D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
2358	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2359	maskvals, namemaskvar, nvm, invardims(4), -1, fractype, &
2360	-1., -1., lairefrac, alaimap)
2361
2362	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_interlai after interpweight_4D'
2363
2364	ENDIF
2365
2366
2367
2368	!
2369	!
2370	!Config Key = RENORM_LAI
2371	!Config Desc = flag to force LAI renormelization
2372	!Config If = LAI_MAP
2373	!Config Def = n
2374	!Config Help = If true, the laimap will be renormalize between llaimin and llaimax parameters.
2375	!Config Units = [FLAG]
2376	!
2377	renormelize_lai = .FALSE.
2378	CALL getin_p('RENORM_LAI',renormelize_lai)
2379
2380	!
2381	laimap(:,:,:) = zero
2382	!
2383	IF (printlev_loc >= 5) THEN
2384	WRITE(numout,*)' slowproc_interlai before starting loop nbpt:', nbpt
2385	END IF
2386
2387	! Assigning the right values and giving a value where information was not found
2388	DO ib=1,nbpt
2389	IF (alaimap(ib) < min_sechiba) THEN
2390	DO jv=1,nvm
2391	laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux
2392	ENDDO
2393	ELSE
2394	DO jv=1, nvm
2395	DO it=1, 12
2396	laimap(ib,jv,it) = lairefrac(ib,jv,it)
2397	ENDDO
2398	ENDDO
2399	END IF
2400	ENDDO
2401	!
2402	! Normelize the read LAI by the values SECHIBA is used to
2403	!
2404	IF ( renormelize_lai ) THEN
2405	DO ib=1,nbpt
2406	DO jv=1, nvm
2407	lmax = MAXVAL(laimap(ib,jv,:))
2408	lmin = MINVAL(laimap(ib,jv,:))
2409	ldelta = lmax-lmin
2410	IF ( ldelta < min_sechiba) THEN
2411	! LAI constante ... keep it constant
2412	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)+(llaimax(jv)+llaimin(jv))/deux
2413	ELSE
2414	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)/(lmax-lmin)*(llaimax(jv)-llaimin(jv))+llaimin(jv)
2415	ENDIF
2416	ENDDO
2417	ENDDO
2418	ENDIF
2419
2420	! Write diagnostics
2421	CALL xios_orchidee_send_field("alaimap",alaimap)
2422	CALL xios_orchidee_send_field("interp_diag_lai",laimap)
2423
2424	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_interlai ended'
2425
2426	END SUBROUTINE slowproc_interlai
2427
2428	!! ================================================================================================================================
2429	!! SUBROUTINE : slowproc_readvegetmax
2430	!!
2431	!>\BRIEF Read and interpolate a vegetation map (by pft)
2432	!!
2433	!! DESCRIPTION : (definitions, functional, design, flags):
2434	!!
2435	!! RECENT CHANGE(S): The subroutine was previously called slowproc_update.
2436	!!
2437	!! MAIN OUTPUT VARIABLE(S):
2438	!!
2439	!! REFERENCE(S) : None
2440	!!
2441	!! FLOWCHART : None
2442	!! \n
2443	!_ ================================================================================================================================
2444
2445	SUBROUTINE slowproc_readvegetmax(nbpt, lalo, neighbours, resolution, contfrac, veget_last, &
2446	veget_next, frac_nobio_next, veget_year, init)
2447
2448	USE interpweight
2449	IMPLICIT NONE
2450
2451	!
2452	!
2453	!
2454	! 0.1 INPUT
2455	!
2456	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs
2457	!! to be interpolated
2458	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
2459	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
2460	!! (1=North and then clockwise)
2461	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
2462	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
2463	!
2464	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last !! old max vegetfrac
2465	INTEGER(i_std), INTENT(in) :: veget_year !! first year for landuse (0 == NO TIME AXIS)
2466	LOGICAL, INTENT(in) :: init !! initialisation : in case of dgvm, it forces update of all PFTs
2467	!
2468	! 0.2 OUTPUT
2469	!
2470	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(out) :: veget_next !! new max vegetfrac
2471	REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next !! new fraction of the mesh which is
2472	!! covered by ice, lakes, ...
2473
2474	!
2475	! 0.3 LOCAL
2476	!
2477	!
2478	CHARACTER(LEN=80) :: filename
2479	INTEGER(i_std) :: ib, inobio, jv
2480	REAL(r_std) :: sumf, err, norm
2481	!
2482	! for DGVM case :
2483	REAL(r_std) :: sum_veg ! sum of vegets
2484	REAL(r_std) :: sum_nobio ! sum of nobios
2485	REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets
2486	REAL(r_std) :: rapport ! (S-B) / (S-A)
2487	LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones)
2488	! e.g. in case of DGVM and not init (optional parameter)
2489	REAL(r_std), DIMENSION(nbpt,nvm) :: vegetrefrac !! veget fractions re-dimensioned
2490	REAL(r_std), DIMENSION(nbpt) :: aveget !! Availability of the soilcol interpolation
2491	REAL(r_std), DIMENSION(nbpt,nvm) :: aveget_nvm !! Availability of the soilcol interpolation
2492	REAL(r_std), DIMENSION(nvm) :: vmin, vmax !! min/max values to use for the renormalization
2493	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2494	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2495	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
2496	!! (variabletypevals(1) = -un, not used)
2497	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2498	!! 'XYKindTime': Input values are kinds
2499	!! of something with a temporal
2500	!! evolution on the dx*dy matrix'
2501	LOGICAL :: nonegative !! whether negative values should be removed
2502	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2503	!! 'nomask': no-mask is applied
2504	!! 'mbelow': take values below maskvals(1)
2505	!! 'mabove': take values above maskvals(1)
2506	!! 'msumrange': take values within 2 ranges;
2507	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2508	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2509	!! (normalized by maskvals(3))
2510	!! 'var': mask values are taken from a
2511	!! variable inside the file (>0)
2512	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2513	!! `maskingtype')
2514	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2515	CHARACTER(LEN=250) :: msg
2516
2517	!_ ================================================================================================================================
2518
2519	IF (printlev_loc >= 5) PRINT *,' In slowproc_readvegetmax'
2520
2521	!
2522	!Config Key = VEGETATION_FILE
2523	!Config Desc = Name of file from which the vegetation map is to be read
2524	!Config If =
2525	!Config Def = PFTmap.nc
2526	!Config Help = The name of the file to be opened to read a vegetation
2527	!Config map (in pft) is to be given here.
2528	!Config Units = [FILE]
2529	!
2530	filename = 'PFTmap.nc'
2531	CALL getin_p('VEGETATION_FILE',filename)
2532	variablename = 'maxvegetfrac'
2533
2534
2535	IF (xios_interpolation) THEN
2536	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readvegetmax: Use XIOS to read and interpolate " &
2537	// TRIM(filename) // " for variable " // TRIM(variablename)
2538
2539	CALL xios_orchidee_recv_field('frac_veget',vegetrefrac)
2540	CALL xios_orchidee_recv_field('frac_veget_frac',aveget_nvm)
2541	aveget(:)=aveget_nvm(:,1)
2542
2543	DO ib = 1, nbpt
2544	IF (aveget(ib) > min_sechiba) THEN
2545	vegetrefrac(ib,:) = vegetrefrac(ib,:)/aveget(ib) ! intersected area normalization
2546	vegetrefrac(ib,:) = vegetrefrac(ib,:)/SUM(vegetrefrac(ib,:))
2547	ENDIF
2548	ENDDO
2549
2550	ELSE
2551
2552	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_readvegetmax: Start interpolate " &
2553	// TRIM(filename) // " for variable " // TRIM(variablename)
2554
2555	! Assigning values to vmin, vmax
2556	vmin = 1
2557	vmax = nvm*1._r_std
2558
2559	variabletypevals = -un
2560
2561	!! Variables for interpweight
2562	! Type of calculation of cell fractions
2563	fractype = 'default'
2564	! Name of the longitude and latitude in the input file
2565	lonname = 'lon'
2566	latname = 'lat'
2567	! Should negative values be set to zero from input file?
2568	nonegative = .FALSE.
2569	! Type of mask to apply to the input data (see header for more details)
2570	maskingtype = 'msumrange'
2571	! Values to use for the masking
2572	maskvals = (/ 1.-1.e-7, 0., 2. /)
2573	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2574	namemaskvar = ''
2575
2576	CALL interpweight_3D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
2577	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2578	maskvals, namemaskvar, nvm, 0, veget_year, fractype, &
2579	-1., -1., vegetrefrac, aveget)
2580	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after interpeeight_3D'
2581	ENDIF
2582	!
2583	! Compute the logical for partial (only anthropic) PTFs update
2584	IF (ok_dgvm .AND. .NOT. init) THEN
2585	partial_update= .TRUE.
2586	ELSE
2587	partial_update=.FALSE.
2588	END IF
2589
2590	IF (printlev_loc >= 5) THEN
2591	WRITE(numout,*)' slowproc_readvegetmax before updating loop nbpt:', nbpt
2592	END IF
2593
2594	IF ( .NOT. partial_update ) THEN
2595	veget_next(:,:)=zero
2596
2597	IF (printlev_loc >=3 .AND. ANY(aveget < min_sechiba)) THEN
2598	WRITE(numout,*) 'Some grid cells on the model grid did not have any points on the source grid.'
2599	IF (init) THEN
2600	WRITE(numout,*) 'Initialization with full fraction of bare soil are done for the below grid cells.'
2601	ELSE
2602	WRITE(numout,*) 'Old values are kept for the below grid cells.'
2603	ENDIF
2604	WRITE(numout,*) 'List of grid cells (ib, lat, lon):'
2605	END IF
2606
2607	DO ib = 1, nbpt
2608	! vegetrefrac is already normalized to sum equal one for each grid cell
2609	veget_next(ib,:) = vegetrefrac(ib,:)
2610
2611	IF (aveget(ib) < min_sechiba) THEN
2612	IF (printlev_loc >=3) WRITE(numout,*) ib,lalo(ib,1),lalo(ib,2)
2613	IF (init) THEN
2614	veget_next(ib,1) = un
2615	veget_next(ib,2:nvm) = zero
2616	ELSE
2617	veget_next(ib,:) = veget_last(ib,:)
2618	ENDIF
2619	ENDIF
2620	ENDDO
2621	ELSE
2622	! Partial update
2623	DO ib = 1, nbpt
2624	IF (aveget(ib) > min_sechiba) THEN
2625	! For the case with properly interpolated grid cells (aveget>0)
2626
2627	! last veget for this point
2628	sum_veg=SUM(veget_last(ib,:))
2629	!
2630	! If the DGVM is activated, only anthropic PFTs are utpdated, the others are copied from previous time-step
2631	veget_next(ib,:) = veget_last(ib,:)
2632
2633	DO jv = 2, nvm
2634	IF ( .NOT. natural(jv) ) THEN
2635	veget_next(ib,jv) = vegetrefrac(ib,jv)
2636	ENDIF
2637	ENDDO
2638
2639	sumvAnthro_old = zero
2640	sumvAnthro = zero
2641	DO jv = 2, nvm
2642	IF ( .NOT. natural(jv) ) THEN
2643	sumvAnthro = sumvAnthro + veget_next(ib,jv)
2644	sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv)
2645	ENDIF
2646	ENDDO
2647
2648	IF ( sumvAnthro_old < sumvAnthro ) THEN
2649	! Increase of non natural vegetations (increase of agriculture)
2650	! The proportion of natural PFT's must be preserved
2651	! ie the sum of vegets is preserved
2652	! and natural PFT / (sum of veget - sum of antropic veget)
2653	! is preserved.
2654	rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old )
2655	DO jv = 1, nvm
2656	IF ( natural(jv) ) THEN
2657	veget_next(ib,jv) = veget_last(ib,jv) * rapport
2658	ENDIF
2659	ENDDO
2660	ELSE
2661	! Increase of natural vegetations (decrease of agriculture)
2662	! The decrease of agriculture is replaced by bare soil. The DGVM will
2663	! re-introduce natural PFT's.
2664	DO jv = 1, nvm
2665	IF ( natural(jv) ) THEN
2666	veget_next(ib,jv) = veget_last(ib,jv)
2667	ENDIF
2668	ENDDO
2669	veget_next(ib,1) = veget_next(ib,1) + sumvAnthro_old - sumvAnthro
2670	ENDIF
2671
2672	! test
2673	IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > 10*EPSILON(un) ) THEN
2674	WRITE(numout,*) 'slowproc_readvegetmax _______'
2675	msg = " No conservation of sum of veget for point "
2676	WRITE(numout,*) TRIM(msg), ib, ",(", lalo(ib,1),",", lalo(ib,2), ")"
2677	WRITE(numout,*) " last sum of veget ", sum_veg, " new sum of veget ", &
2678	SUM(veget_next(ib,:)), " error : ", SUM(veget_next(ib,:))-sum_veg
2679	WRITE(numout,*) " Anthropic modifications : last ",sumvAnthro_old," new ",sumvAnthro
2680	CALL ipslerr_p (3,'slowproc_readvegetmax', &
2681	& 'No conservation of sum of veget_next', &
2682	& "The sum of veget_next is different after reading Land Use map.", &
2683	& '(verify the dgvm case model.)')
2684	ENDIF
2685	ELSE
2686	! For the case when there was a propblem with the interpolation, aveget < min_sechiba
2687	WRITE(numout,*) 'slowproc_readvegetmax _______'
2688	WRITE(numout,*) " No land point in the map for point ", ib, ",(", lalo(ib,1), ",", &
2689	lalo(ib,2),")"
2690	CALL ipslerr_p (2,'slowproc_readvegetmax', &
2691	& 'Problem with vegetation file for Land Use.', &
2692	& "No land point in the map for point", &
2693	& '(verify your land use file.)')
2694	veget_next(ib,:) = veget_last(ib,:)
2695	ENDIF
2696
2697	ENDDO
2698	ENDIF
2699
2700	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after updating'
2701	!
2702	frac_nobio_next (:,:) = un
2703	!
2704	!MM
2705	! Work only for one nnobio !! (ie ice)
2706	DO inobio=1,nnobio
2707	DO jv=1,nvm
2708	DO ib = 1, nbpt
2709	frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv)
2710	ENDDO
2711	ENDDO
2712	ENDDO
2713
2714	DO ib = 1, nbpt
2715	sum_veg = SUM(veget_next(ib,:))
2716	sum_nobio = SUM(frac_nobio_next(ib,:))
2717	IF (sum_nobio < 0.) THEN
2718	frac_nobio_next(ib,:) = zero
2719	veget_next(ib,1) = veget_next(ib,1) + sum_nobio
2720	sum_veg = SUM(veget_next(ib,:))
2721	ENDIF
2722	sumf = sum_veg + sum_nobio
2723	IF (sumf > min_sechiba) THEN
2724	veget_next(ib,:) = veget_next(ib,:) / sumf
2725	frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf
2726	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2727	err=norm-un
2728	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ",ib, &
2729	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un, sumf",err,sumf
2730	IF (abs(err) > -EPSILON(un)) THEN
2731	IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN
2732	frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err
2733	ELSE
2734	veget_next(ib,1) = veget_next(ib,1) - err
2735	ENDIF
2736	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2737	err=norm-un
2738	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ", ib, &
2739	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un",err
2740	IF (abs(err) > EPSILON(un)) THEN
2741	WRITE(numout,*) ' slowproc_readvegetmax _______'
2742	WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2743	WRITE(numout,*) " err(sum-1.) = ",abs(err)
2744	CALL ipslerr_p (2,'slowproc_readvegetmax', &
2745	& 'Problem with sum vegetation + sum fracnobio for Land Use.', &
2746	& "sum not equal to 1.", &
2747	& '(verify your land use file.)')
2748	aveget(ib) = -0.6
2749	ENDIF
2750	ENDIF
2751	ELSE
2752	! sumf < min_sechiba
2753	WRITE(numout,*) ' slowproc_readvegetmax _______'
2754	WRITE(numout,*)" No vegetation nor frac_nobio for point ", ib, ",(", lalo(ib,1), ",", &
2755	lalo(ib,2),")"
2756	WRITE(numout,*)" Replaced by bare_soil !! "
2757	veget_next(ib,1) = un
2758	veget_next(ib,2:nvm) = zero
2759	frac_nobio_next(ib,:) = zero
2760	!!!$ CALL ipslerr_p (3,'slowproc_readvegetmax', &
2761	!!!$ & 'Problem with vegetation file for Land Use.', &
2762	!!!$ & "No vegetation nor frac_nobio for point ", &
2763	!!!$ & '(verify your land use file.)')
2764	ENDIF
2765	ENDDO
2766
2767	!! Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
2768	!! Normalize to have the sum equal 1.
2769	CALL slowproc_veget_max_limit(nbpt, frac_nobio_next, veget_next)
2770
2771	! Write diagnostics
2772	CALL xios_orchidee_send_field("aveget",aveget)
2773	CALL xios_orchidee_send_field("interp_diag_aveget",aveget)
2774	CALL xios_orchidee_send_field("interp_diag_vegetrefrac",vegetrefrac)
2775
2776	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_readvegetmax ended'
2777
2778	END SUBROUTINE slowproc_readvegetmax
2779
2780
2781	!! ================================================================================================================================
2782	!! SUBROUTINE : slowproc_nearest
2783	!!
2784	!>\BRIEF looks for nearest grid point on the fine map
2785	!!
2786	!! DESCRIPTION : (definitions, functional, design, flags):
2787	!!
2788	!! RECENT CHANGE(S): None
2789	!!
2790	!! MAIN OUTPUT VARIABLE(S): ::inear
2791	!!
2792	!! REFERENCE(S) : None
2793	!!
2794	!! FLOWCHART : None
2795	!! \n
2796	!_ ================================================================================================================================
2797
2798	SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
2799
2800	!! INTERFACE DESCRIPTION
2801
2802	!! 0.1 input variables
2803
2804	INTEGER(i_std), INTENT(in) :: iml !! size of the vector
2805	REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5 !! longitude and latitude vector, for the 5km vegmap
2806	REAL(r_std), INTENT(in) :: lonmod, latmod !! longitude and latitude modelled
2807
2808	!! 0.2 output variables
2809
2810	INTEGER(i_std), INTENT(out) :: inear !! location of the grid point from the 5km vegmap grid
2811	!! closest from the modelled grid point
2812
2813	!! 0.4 Local variables
2814
2815	REAL(r_std) :: pa, p
2816	REAL(r_std) :: coscolat, sincolat
2817	REAL(r_std) :: cospa, sinpa
2818	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
2819	INTEGER(i_std) :: i
2820	INTEGER(i_std), DIMENSION(1) :: ineartab
2821	INTEGER :: ALLOC_ERR
2822
2823	!_ ================================================================================================================================
2824
2825	ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
2826	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_nearest','Error in allocation for cosang','','')
2827
2828	pa = pi/2.0 - latmod*pi/180.0 ! dist. between north pole and the point a
2829	!! COLATITUDE, in radian
2830	cospa = COS(pa)
2831	sinpa = SIN(pa)
2832
2833	DO i = 1, iml
2834
2835	sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 ) !! sinus of the colatitude
2836	coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 ) !! cosinus of the colatitude
2837
2838	p = (lonmod-lon5(i))*pi/180.0 !! angle between a & b (between their meridian)in radians
2839
2840	!! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
2841	cosang(i) = cospacoscolat + sinpasincolat*COS(p) !! TL : cosang is maximum when angle is at minimal value
2842	!! orthodromic distance between 2 points : cosang = cosinus (arc(AB)/R), with
2843	!R = Earth radius, then max(cosang) = max(cos(arc(AB)/R)), reached when arc(AB)/R is minimal, when
2844	! arc(AB) is minimal, thus when point B (corresponding grid point from LAI MAP) is the nearest from
2845	! modelled A point
2846	ENDDO
2847
2848	ineartab = MAXLOC( cosang(:) )
2849	inear = ineartab(1)
2850
2851	DEALLOCATE(cosang)
2852	END SUBROUTINE slowproc_nearest
2853
2854	!! ================================================================================================================================
2855	!! SUBROUTINE : slowproc_soilt
2856	!!
2857	!>\BRIEF Interpolate the Zobler or Reynolds/USDA soil type map
2858	!!
2859	!! DESCRIPTION : (definitions, functional, design, flags):
2860	!!
2861	!! RECENT CHANGE(S): Nov 2014, ADucharne
2862	!! Nov 2020, Salma Tafasca and Agnes Ducharne: adding a choice for spmipexp/SPMIPEXP,
2863	!! and everything needed to read all maps and assign parameter values.
2864	!!
2865	!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction
2866	!!
2867	!! REFERENCE(S) : Reynold, Jackson, and Rawls (2000). Estimating soil water-holding capacities
2868	!! by linking the Food and Agriculture Organization soil map of the world with global pedon
2869	!! databases and continuous pedotransfer functions, WRR, 36, 3653-3662
2870	!!
2871	!! FLOWCHART : None
2872	!! \n
2873	!_ ================================================================================================================================
2874	SUBROUTINE slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction)
2875
2876	USE interpweight
2877
2878	IMPLICIT NONE
2879	!
2880	!
2881	! This subroutine should read the Zobler/Reynolds map and interpolate to the model grid.
2882	! The method is to get fraction of the three/12 main soiltypes for each grid box.
2883	! For the Zobler case, also called FAO in the code, the soil fraction are going to be put
2884	! into the array soiltype in the following order : coarse, medium and fine.
2885	! For the Reynolds/USDA case, the soiltype array follows the order defined in constantes_soil_var.f90
2886	!
2887	!
2888	!! 0.1 INPUT
2889	!
2890	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
2891	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
2892	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2893	!! (1=North and then clockwise)
2894	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
2895	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
2896	!
2897	! 0.2 OUTPUT
2898	!
2899	!salma: added soil params and njsc because needed in the calculation of the soil params
2900	INTEGER(i_std),DIMENSION (nbpt), INTENT (out) :: njsc !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
2901	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: ks !! Hydraulic conductivity at saturation (mm {-1})
2902	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: nvan !! Van Genuchten coeficients n (unitless)
2903	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: avan !! Van Genuchten coeficients a (mm-1})
2904	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: mcr !! Residual volumetric water content (m^{3} m^{-3})
2905	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: mcs !! Saturated volumetric water content (m^{3} m^{-3})
2906	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: mcfc !! Volumetric water content at field capacity (m^{3} m^{-3})
2907	REAL(r_std),DIMENSION (nbpt), INTENT (out) :: mcw !! Volumetric water content at wilting point (m^{3} m^{-3})
2908
2909	REAL(r_std), INTENT(out) :: soilclass(nbpt, nscm) !! Soil type map to be created from the Zobler map
2910	!! or a map defining the 12 USDA classes (e.g. Reynolds)
2911	!! Holds the area of each texture class in the ORCHIDEE grid cells
2912	!! Final unit = fraction of ORCHIDEE grid-cell (unitless)
2913	REAL(r_std), INTENT(out) :: clayfraction(nbpt) !! The fraction of clay as used by STOMATE
2914	REAL(r_std), INTENT(out) :: sandfraction(nbpt) !! The fraction of sand (for SP-MIP)
2915	REAL(r_std), INTENT(out) :: siltfraction(nbpt) !! The fraction of silt (for SP-MIP)
2916	!
2917	!
2918	! 0.3 LOCAL
2919	!
2920	!salma: added the following local variable to be used for all the soil hydraulic parameters
2921	REAL(r_std), DIMENSION(nbpt) :: param !! to be introduced in function: interpweight
2922
2923	CHARACTER(LEN=80) :: filename
2924	INTEGER(i_std) :: ib, ilf, nbexp, i
2925	INTEGER(i_std) :: fopt !! Nb of pts from the texture map within one ORCHIDEE grid-cell
2926	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt !! Texture the different points from the input texture map
2927	!! in one ORCHIDEE grid cell (unitless)
2928	!
2929	! Number of texture classes in Zobler
2930	!
2931	INTEGER(i_std), PARAMETER :: nzobler = 7 !! Nb of texture classes according in the Zobler map
2932	REAL(r_std),ALLOCATABLE :: textfrac_table(:,:) !! conversion table between the texture index
2933	!! and the granulometric composition
2934	!
2935	INTEGER :: ALLOC_ERR
2936	INTEGER :: ntextinfile !! number of soil textures in the in the file
2937	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: textrefrac !! text fractions re-dimensioned
2938	REAL(r_std), DIMENSION(nbpt) :: atext !! Availability of the texture interpolation
2939	!salma added the following 3 variables to control the SP-MIP simulations
2940	REAL(r_std), DIMENSION(nbpt) :: aparam !! Availability of the parameter interpolation
2941	CHARACTER(LEN=80) :: spmipexp !! designing the number of sp-mip experiment
2942	CHARACTER(LEN=80) :: unif_case !! designing the model of experiment 4 (sp_mip)
2943
2944	REAL(r_std) :: vmin, vmax !! min/max values to use for the
2945
2946	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2947	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in input file
2948	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
2949	!! (variabletypevals(1) = -un, not used)
2950	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2951	!! 'XYKindTime': Input values are kinds
2952	!! of something with a temporal
2953	!! evolution on the dx*dy matrix'
2954	LOGICAL :: nonegative !! whether negative values should be removed
2955	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2956	!! 'nomask': no-mask is applied
2957	!! 'mbelow': take values below maskvals(1)
2958	!! 'mabove': take values above maskvals(1)
2959	!! 'msumrange': take values within 2 ranges;
2960	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2961	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2962	!! (normalized by maskvals(3))
2963	!! 'var': mask values are taken from a
2964	!! variable inside the file (>0)
2965	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2966	!! `maskingtype')
2967	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2968	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
2969	REAL(r_std) :: sgn !! sum of fractions excluding glaciers and ocean
2970
2971	! For the calculation of field capacity and wilting point
2972	REAL(r_std),DIMENSION (nbpt) :: mvan !! Van Genuchten parameter m
2973	REAL(r_std),DIMENSION (nbpt) :: psi_w !! Matrix potential characterizing the wilting point (mm)
2974	REAL(r_std),DIMENSION (nbpt) :: psi_fc !! Matrix potential characterizing the field capacity (mm)
2975
2976	!_ ================================================================================================================================
2977
2978	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt'
2979
2980	! The soil parameters are defined by several keywords in run.def:
2981	! (a) soil_classif tells which kind of soil texture map you will read (mandatory):
2982	! - usda for 12 USDA texture classes (Reynolds, SoilGrids, SPMIP, etc) updated to 13 classes
2983	! for clay oxisols by Salma Tafasca
2984	! - zobler to read teh Zobler map and reduce it to 3 classes (fine, medium, coarse)
2985	! (b) spmipexp was introduced by Salma Tafasca for the SPMIP project
2986	! maps: Reading the soil parameter maps of SPMIP
2987	! unif: Imposing uniform soil texture over the globe (4 texture options, with parameter values imposed by SP-MIP)
2988	! Even with maps, some parameters (thermics) are defined based on texture.
2989	! So we read a soil texture map in all experiments but unif, where soil texture is imposed by njsc(:).
2990	! (c) unif_case to choose the soil texture assigned if spmipexp=maps (4 hard_coded possibilities)
2991
2992	! IMPORTANT: if no spmipexp is defined in run.def, the model works as before, by deriving the soil parameters
2993	! from a soil texture map, itself defined by the SOILTYPE_CLASSIF keyword, and soil_classif variable
2994	! 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
2995
2996	!Config Key = SPMIPEXP
2997	!Config Desc = Types of alternative hydraulic parameters
2998	!Config Def = 'texture'
2999	!Config If =
3000	!Config Help = possible values: maps, unif
3001	!Config Units = [-]
3002	spmipexp='texture' ! default is to define parameters from soil texture, with soil_classif = 'zobler' or 'usda'
3003	CALL getin_p("SPMIPEXP",spmipexp)
3004
3005	IF (spmipexp == 'unif') THEN
3006	! case where unif=exp4 is selected: uniform soil parameters
3007	! the values of the hydraulic parameters below come from SP-MIP,
3008	! and correspond to the Rosetta PTF (Schaap et al., 2001)
3009
3010	! sp_mip_experiment_4: select another level of experiment: a, b, c or d in run.def
3011
3012	!Config Key = UNIF_CASE
3013	!Config Desc = Types of uniform soil textures in SPMIP
3014	!Config Def = 'b'
3015	!Config If = SPMIPEXP='unif'
3016	!Config Help = possible values: a, b, c and d
3017	!Config Units = [-]
3018	unif_case='b' ! default = loamy soil
3019	CALL getin_p("UNIF_CASE",unif_case)
3020
3021	SELECTCASE (unif_case)
3022
3023	CASE ('a') ! loamy sand
3024	clayfraction=0.06
3025	sandfraction=0.81
3026	siltfraction=0.13
3027	DO ib=1 , nbpt
3028	njsc(ib) = 2
3029	mcr(ib) = 0.049
3030	mcs(ib) = 0.39
3031	ks(ib) = (1.41e-5)100024*3600
3032	avan(ib) = 3.475*(1e-3)
3033	nvan(ib) = 1.746
3034	mcfc(ib) = 0.1039
3035	mcw(ib) = 0.05221
3036	ENDDO
3037
3038	CASE ('b') !loam
3039	clayfraction=0.2
3040	sandfraction=0.4
3041	siltfraction=0.4
3042	DO ib=1, nbpt
3043	njsc(ib) = 6
3044	mcr(ib) = 0.061
3045	mcs(ib) = 0.399
3046	ks(ib) = (3.38e-6)100024*3600
3047	avan(ib) = 1.112*(1e-3)
3048	nvan(ib) = 1.472
3049	mcfc(ib) = 0.236
3050	mcw(ib) = 0.09115
3051	ENDDO
3052
3053	CASE ('c') !silt
3054	clayfraction=0.1
3055	sandfraction=0.06
3056	siltfraction=0.84
3057	DO ib=1, nbpt
3058	njsc(ib)=5
3059	mcr(ib) = 0.05
3060	mcs(ib) = 0.489
3061	ks(ib) = (2.81e-6)100024*3600
3062	avan(ib) = 0.6577*(1e-3)
3063	nvan(ib) = 1.679
3064	mcfc(ib) = 0.2854
3065	mcw(ib) = 0.06944
3066	ENDDO
3067
3068	CASE ('d')!clay
3069	clayfraction=0.55
3070	sandfraction=0.15
3071	siltfraction=0.3
3072	DO ib=1, nbpt
3073	njsc(ib)=12
3074	mcr(ib) = 0.098
3075	mcs(ib) = 0.459
3076	ks(ib) = (9.74e-7)100024*3600
3077	avan(ib) = 1.496*(1e-3)
3078	nvan(ib) = 1.253
3079	mcfc(ib) = 0.3329
3080	mcw(ib) = 0.1897
3081	ENDDO
3082
3083	CASE DEFAULT
3084
3085	WRITE (numout,*) 'Unsupported experiment number. Choose between a, b, c or d according to sp_mip_experiment_4 number'
3086	CALL ipslerr_p(3,'hydrol_init','Unsupported experiment number. ',&
3087	'Choose between a,b,c or d','')
3088	ENDSELECT
3089
3090	ELSE ! spmipexp is either exp1=maps, or texture for exp2 or exp3 (or typing error!)
3091
3092	! In these cases (maps or texture), we need to read the soil texture map
3093
3094	!Config Key = SOILCLASS_FILE
3095	!Config Desc = Name of file from which soil types are read
3096	!Config Def = soils_param.nc
3097	!Config If = NOT(IMPOSE_VEG)
3098	!Config Help = The name of the file to be opened to read the soil types.
3099	!Config The data from this file is then interpolated to the grid of
3100	!Config of the model. The aim is to get fractions for sand loam and
3101	!Config clay in each grid box. This information is used for soil hydrology
3102	!Config and respiration.
3103	!Config Units = [FILE]
3104	!
3105	! soils_param.nc file is 1deg soil texture file (Zobler)
3106	! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution)
3107
3108	filename = 'soils_param.nc'
3109	CALL getin_p('SOILCLASS_FILE',filename)
3110
3111	variablename = 'soiltext'
3112
3113	!! Variables for interpweight
3114	! Type of calculation of cell fractions
3115	fractype = 'default'
3116
3117	IF (xios_interpolation) THEN
3118	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " &
3119	// TRIM(filename) // " for variable " // TRIM(variablename)
3120
3121	SELECT CASE(soil_classif)
3122
3123	CASE('none')
3124	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
3125	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3126	DO ib=1, nbpt
3127	njsc(ib) = usda_default ! 6 = Loam
3128	clayfraction(ib) = clayfrac_usda(usda_default)
3129	sandfraction(ib) = sandfrac_usda(usda_default)
3130	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3131	ENDDO
3132
3133	CASE('zobler')
3134	! !
3135	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification, to be read using XIOS"
3136	!
3137	ALLOCATE(textrefrac(nbpt,nzobler))
3138	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
3139	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3140	CALL get_soilcorr_zobler (nzobler, textfrac_table)
3141	CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
3142	CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
3143	CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
3144	CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
3145	CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
3146	CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
3147	CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
3148
3149	CALL get_soilcorr_zobler (nzobler, textfrac_table)
3150	! !
3151	DO ib =1, nbpt
3152	soilclass(ib,:)=0.
3153	soilclass(ib,fao2usda(1))=textrefrac(ib,1)
3154	soilclass(ib,fao2usda(2))=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
3155	soilclass(ib,fao2usda(3))=textrefrac(ib,5)
3156
3157	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
3158	! over the zobler pixels composing the ORCHIDEE grid-cell
3159	clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + &
3160	textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + &
3161	textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7)
3162
3163	sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + &
3164	textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + &
3165	textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7)
3166
3167	siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + &
3168	textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + &
3169	textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7)
3170
3171	sgn=SUM(soilclass(ib,:)) ! grid-cell fraction with texture info
3172
3173	IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume that texture = Loam
3174	njsc(ib) = usda_default ! 6 = Loam
3175	clayfraction(ib) = clayfrac_usda(usda_default)
3176	sandfraction(ib) = sandfrac_usda(usda_default)
3177	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3178	atext(ib)=0.
3179	ELSE
3180	atext(ib)=sgn
3181	clayfraction(ib) = clayfraction(ib) / sgn
3182	sandfraction(ib) = sandfraction(ib) / sgn
3183	siltfraction(ib) = siltfraction(ib) / sgn
3184	soilclass(ib,:) = soilclass(ib,:) / sgn
3185	njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class
3186	ENDIF
3187
3188	ENDDO
3189
3190	CASE('usda')
3191
3192	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda'
3193	!
3194	WRITE(numout,*) "Using a soilclass map with usda classification, to be read using XIOS"
3195	!
3196	ALLOCATE(textrefrac(nbpt,nscm))
3197	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
3198	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3199
3200	CALL get_soilcorr_usda (nscm, textfrac_table)
3201
3202	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
3203
3204	CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
3205	CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
3206	CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
3207	CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
3208	CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
3209	CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
3210	CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
3211	CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8))
3212	CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9))
3213	CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10))
3214	CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11))
3215	CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12))
3216	CALL xios_orchidee_recv_field('soiltext13',textrefrac(:,13))
3217
3218	CALL get_soilcorr_usda (nscm, textfrac_table)
3219	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
3220
3221	DO ib =1, nbpt
3222	clayfraction(ib) = 0.0
3223	DO ilf = 1,nscm
3224	soilclass(ib,ilf)=textrefrac(ib,ilf)
3225	clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf)
3226	sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf)
3227	siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf)
3228	! textfrac_table holds the %silt,%sand,%clay
3229	ENDDO
3230
3231	sgn=SUM(soilclass(ib,:)) ! grid-cell fraction with texture info
3232
3233	IF (sgn < min_sechiba) THEN ! if no texture info in this grid-point, we assume that texture = Loam
3234	njsc(ib) = usda_default ! 6 = Loam
3235	clayfraction(ib) = clayfrac_usda(usda_default)
3236	sandfraction(ib) = sandfrac_usda(usda_default)
3237	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3238	atext(ib)=0
3239	ELSE
3240	soilclass(ib,:) = soilclass(ib,:) / sgn
3241	clayfraction(ib) = clayfraction(ib) / sgn
3242	sandfraction(ib) = sandfraction(ib) / sgn
3243	siltfraction(ib) = siltfraction(ib) / sgn
3244	atext(ib)=sgn
3245	njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class
3246	ENDIF
3247
3248	ENDDO
3249
3250	CASE DEFAULT
3251	WRITE(numout,*) 'slowproc_soilt:'
3252	WRITE(numout,*) ' A non supported soil type classification has been chosen'
3253	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
3254	END SELECT
3255
3256	ELSE ! xios_interpolation
3257	! Read and interpolate using stardard method with IOIPSL and aggregate
3258
3259	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
3260	// TRIM(filename) // " for variable " // TRIM(variablename)
3261
3262	! Name of the longitude and latitude in the input file
3263	lonname = 'nav_lon'
3264	latname = 'nav_lat'
3265
3266	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " &
3267	// TRIM(filename) // " for variable " // TRIM(variablename)
3268
3269	IF ( TRIM(soil_classif) /= 'none' ) THEN
3270
3271	! Define a variable for the number of soil textures in the input file
3272	SELECTCASE(soil_classif)
3273	CASE('zobler')
3274	ntextinfile=nzobler
3275	CASE('usda')
3276	ntextinfile=nscm
3277	CASE DEFAULT
3278	WRITE(numout,*) 'slowproc_soilt:'
3279	WRITE(numout,*) ' A non supported soil type classification has been chosen'
3280	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
3281	ENDSELECT
3282
3283	ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR)
3284	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',&
3285	'','')
3286
3287	! Assigning values to vmin, vmax
3288	vmin = un
3289	vmax = ntextinfile*un
3290
3291	ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR)
3292	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','')
3293	variabletypevals = -un
3294
3295	!! Variables for interpweight
3296	! Should negative values be set to zero from input file?
3297	nonegative = .FALSE.
3298	! Type of mask to apply to the input data (see header for more details)
3299	maskingtype = 'mabove'
3300	! Values to use for the masking
3301	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3302	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
3303
3304	namemaskvar = ''
3305
3306	CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours, &
3307	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3308	maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext)
3309
3310	ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR)
3311	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','')
3312	ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR)
3313	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','')
3314
3315	IF (printlev_loc >= 5) THEN
3316	WRITE(numout,*)' slowproc_soilt after interpweight_2D'
3317	WRITE(numout,*)' slowproc_soilt before starting loop nbpt:', nbpt
3318	WRITE(numout,*)" slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..."
3319	END IF
3320	ELSE
3321	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt using default values all points are propertly ' // &
3322	'interpolated atext = 1. everywhere!'
3323	atext = 1.
3324	END IF
3325
3326	nbexp = 0
3327	SELECTCASE(soil_classif)
3328	CASE('none')
3329	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
3330	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3331	DO ib=1, nbpt
3332	njsc(ib) = usda_default ! 6 = Loam
3333	clayfraction(ib) = clayfrac_usda(usda_default)
3334	sandfraction(ib) = sandfrac_usda(usda_default)
3335	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3336	ENDDO
3337	CASE('zobler')
3338	! !
3339	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
3340	!
3341	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
3342	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3343	CALL get_soilcorr_zobler (nzobler, textfrac_table)
3344
3345	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt after getting table of textures'
3346	DO ib =1, nbpt
3347	soilclass(ib,:) = zero
3348	clayfraction(ib) = zero
3349	sandfraction(ib) = zero
3350	siltfraction(ib) = zero
3351	!
3352	! vecpos: List of positions where textures were not zero
3353	! vecpos(1): number of not null textures found
3354	vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq')
3355	fopt = vecpos(1)
3356
3357	IF ( fopt .EQ. 0 ) THEN
3358	! No points were found for current grid box, use default values
3359	nbexp = nbexp + 1
3360	njsc(ib) = usda_default ! 6=Loam
3361	clayfraction(ib) = clayfrac_usda(usda_default)
3362	sandfraction(ib) = sandfrac_usda(usda_default)
3363	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3364
3365	ELSE
3366	IF (fopt == nzobler) THEN
3367	! All textures are not zero
3368	solt=(/(i,i=1,nzobler)/)
3369	ELSE
3370	DO ilf = 1,fopt
3371	solt(ilf) = vecpos(ilf+1)
3372	END DO
3373	END IF
3374	!
3375	! Compute the fraction of each textural class
3376	!
3377	sgn = 0.
3378	DO ilf = 1,fopt
3379	!
3380	! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE
3381	! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine)
3382	! For type 6 = glacier, default values are set and it is also taken into account during the normalization
3383	! of the fractions (done in interpweight_2D)
3384	! Note that type 0 corresponds to ocean but it is already removed using the mask above.
3385	!
3386	IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. &
3387	(solt(ilf) .NE. 6) ) THEN
3388	SELECT CASE(solt(ilf))
3389	CASE(1)
3390	soilclass(ib,fao2usda(1)) = soilclass(ib,fao2usda(1)) + textrefrac(ib,solt(ilf))
3391	CASE(2)
3392	soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
3393	CASE(3)
3394	soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
3395	CASE(4)
3396	soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
3397	CASE(5)
3398	soilclass(ib,fao2usda(3)) = soilclass(ib,fao2usda(3)) + textrefrac(ib,solt(ilf))
3399	CASE(7)
3400	soilclass(ib,fao2usda(2)) = soilclass(ib,fao2usda(2)) + textrefrac(ib,solt(ilf))
3401	CASE DEFAULT
3402	WRITE(numout,*) 'We should not be here, an impossible case appeared'
3403	CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','')
3404	END SELECT
3405	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
3406	! over the zobler pixels composing the ORCHIDEE grid-cell
3407	clayfraction(ib) = clayfraction(ib) + &
3408	& textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf))
3409	sandfraction(ib) = sandfraction(ib) + &
3410	& textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf))
3411	siltfraction(ib) = siltfraction(ib) + &
3412	& textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf))
3413	! Sum the fractions which are not glaciers nor ocean
3414	sgn = sgn + textrefrac(ib,solt(ilf))
3415	ELSE
3416	IF (solt(ilf) .GT. nzobler) THEN
3417	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
3418	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','')
3419	ENDIF
3420	END IF
3421	ENDDO
3422
3423	IF ( sgn .LT. min_sechiba) THEN
3424	! Set default values if grid cells were only covered by glaciers or ocean
3425	! or if now information on the source grid was found.
3426	nbexp = nbexp + 1
3427	njsc(ib) = usda_default ! 6 = Loam
3428	clayfraction(ib) = clayfrac_usda(usda_default)
3429	sandfraction(ib) = sandfrac_usda(usda_default)
3430	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3431	ELSE
3432	! Normalize using the fraction of surface not including glaciers and ocean
3433	soilclass(ib,:) = soilclass(ib,:)/sgn
3434	clayfraction(ib) = clayfraction(ib)/sgn
3435	sandfraction(ib) = sandfraction(ib)/sgn
3436	siltfraction(ib) = siltfraction(ib)/sgn
3437	njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class
3438	ENDIF
3439	ENDIF
3440	ENDDO
3441
3442	! The "USDA" case reads a map of the 12 USDA texture classes,
3443	! such as to assign the corresponding soil properties
3444	CASE("usda")
3445	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification"
3446
3447	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
3448	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3449
3450	CALL get_soilcorr_usda (nscm, textfrac_table)
3451
3452	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
3453	!
3454	DO ib =1, nbpt
3455	! GO through the point we have found
3456	!
3457	! Provide which textures were found
3458	! vecpos: List of positions where textures were not zero
3459	! vecpos(1): number of not null textures found
3460	vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq')
3461	fopt = vecpos(1)
3462	!
3463	! Check that we found some points
3464	!
3465	soilclass(ib,:) = 0.0
3466	clayfraction(ib) = 0.0
3467	sandfraction(ib) = 0.0
3468	siltfraction(ib) = 0.0
3469
3470	IF ( fopt .EQ. 0) THEN
3471	! No points were found for current grid box, use default values
3472	IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib
3473	nbexp = nbexp + 1
3474	njsc(ib) = usda_default ! 6 = Loam
3475	clayfraction(ib) = clayfrac_usda(usda_default)
3476	sandfraction(ib) = sandfrac_usda(usda_default)
3477	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3478	ELSE
3479	IF (fopt == nscm) THEN
3480	! All textures are not zero
3481	solt(:) = (/(i,i=1,nscm)/)
3482	ELSE
3483	DO ilf = 1,fopt
3484	solt(ilf) = vecpos(ilf+1)
3485	END DO
3486	END IF
3487
3488	! Compute the fraction of each textural class
3489	DO ilf = 1,fopt
3490	IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
3491	soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf))
3492	clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * &
3493	textrefrac(ib,solt(ilf))
3494	sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * &
3495	textrefrac(ib,solt(ilf))
3496	siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * &
3497	textrefrac(ib,solt(ilf))
3498	ELSE
3499	IF (solt(ilf) .GT. nscm) THEN
3500	WRITE(,) 'The file contains a soil color class which is incompatible with this program'
3501	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','')
3502	ENDIF
3503	ENDIF
3504	!
3505	ENDDO
3506	njsc(ib) = MAXLOC(soilclass(ib,:),1) ! Dominant texture class
3507
3508	! Set default values if the surface in source file is too small
3509	! Warning - This test is donne differently for Zobler (based on sgn, related to class 6=ice)
3510	IF ( atext(ib) .LT. min_sechiba) THEN
3511	nbexp = nbexp + 1
3512	njsc(ib) = usda_default ! 6 = Loam
3513	clayfraction(ib) = clayfrac_usda(usda_default)
3514	sandfraction(ib) = sandfrac_usda(usda_default)
3515	siltfraction(ib) = 1.-clayfrac_usda(usda_default)-sandfrac_usda(usda_default)
3516	ENDIF
3517	ENDIF
3518
3519	ENDDO
3520
3521	IF (printlev_loc>=4) WRITE (numout,*) ' slowproc_soilt: End case usda'
3522
3523	CASE DEFAULT
3524	WRITE(numout,*) 'slowproc_soilt _______'
3525	WRITE(numout,*) ' A non supported soil type classification has been chosen'
3526	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
3527	ENDSELECT
3528	IF (printlev_loc >= 5 ) WRITE(numout,*)' slowproc_soilt end of type classification'
3529
3530	IF ( nbexp .GT. 0 ) THEN
3531	WRITE(numout,*) 'slowproc_soilt:'
3532	WRITE(numout,*) ' The interpolation of variable soiltext had ', nbexp
3533	WRITE(numout,*) ' points without data. This are either coastal points or ice covered land.'
3534	WRITE(numout,*) ' The problem was solved by using the default soil types.'
3535	ENDIF
3536
3537	IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals)
3538	IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac)
3539	IF (ALLOCATED(solt)) DEALLOCATE (solt)
3540	IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table)
3541
3542	ENDIF ! xios_interpolation
3543
3544	! End of soil texture reading, for 'maps' and classical behavior
3545
3546	IF (spmipexp == 'maps') THEN
3547	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt: Read soil hydraulic parameters with IOIPSL'
3548
3549	! Read using IOIPSL and interpolate using aggregate tool in ORCHIDEE
3550
3551	!Config Key = PARAM_FILE
3552	!Config Desc = Name of file from which soil parameter values are read
3553	!Config Def = params_sp_mip.nc
3554	!Config If = params_sp_mip.nc
3555	!Config Help = The name of the file to be opened to read values of parameters.
3556	!Config The data from this file is then interpolated to the grid of
3557	!Config of the model.
3558	!Config Units = [FILE]
3559	!
3560	! params_sp_mip.nc file is 0.5 deg soil hydraulic parameters file provided by sp_mip
3561
3562	filename = 'params_sp_mip.nc'
3563	CALL getin_p('PARAM_FILE',filename)
3564
3565	!! Variables for interpweight
3566	! Type of calculation of cell fractions
3567	fractype = 'default'
3568	! Name of the longitude and latitude in the input file
3569	lonname = 'nav_lon'
3570	latname = 'nav_lat'
3571	! Assigning values to vmin, vmax (there are not types/categories
3572	vmin =0.
3573	vmax = 99999.
3574	!! Variables for interpweight
3575	! Should negative values be set to zero from input file?
3576	nonegative = .FALSE.
3577	! Type of mask to apply to the input data (see header for more details)
3578	maskingtype = 'mabove'
3579	! Values to use for the masking
3580	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3581	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
3582	namemaskvar = ''
3583
3584	variablename = 'ks'
3585	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
3586	// TRIM(filename) // " for variable " // TRIM(variablename)
3587	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3588	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3589	maskvals, namemaskvar, -1, fractype, 0., 0., &
3590	ks, aparam)
3591	WRITE(numout,*) 'ks map is read _______'
3592
3593	variablename = 'alpha'
3594	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3595	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3596	maskvals, namemaskvar, -1, fractype, 0., 0., &
3597	avan, aparam)
3598	WRITE(numout,*) 'avan map read _______'
3599
3600	variablename = 'thetar'
3601	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3602	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3603	maskvals, namemaskvar, -1, fractype, 0., 0., &
3604	mcr, aparam)
3605	WRITE(numout,*) 'thetar map read _______'
3606
3607	variablename = 'thetas'
3608	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3609	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3610	maskvals, namemaskvar, -1, fractype, 0., 0., &
3611	mcs, aparam)
3612	WRITE(numout,*) 'thetas map read _______'
3613
3614	variablename = 'thetapwpvg' ! mcw
3615	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3616	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3617	maskvals, namemaskvar, -1, fractype, 0., 0., &
3618	mcw, aparam)
3619	WRITE(numout,*) 'thetapwpvg map read _______'
3620
3621	variablename = 'thetafcvg' !mcfc
3622	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3623	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3624	maskvals, namemaskvar, -1, fractype, 0., 0., &
3625	mcfc, aparam)
3626	WRITE(numout,*) 'thetafcvg map read _______'
3627
3628	variablename = 'nvg'
3629	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3630	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3631	maskvals, namemaskvar, -1, fractype, 0., 0., &
3632	nvan, aparam)
3633	WRITE(numout,*) 'nvan map read _______'
3634
3635	ELSE ! spmipexp is not maps nor unif, then it must be texture
3636	IF (spmipexp == 'texture') THEN
3637	! Whichever the soil texture map, we can use the USDA parameter vectors with 13 values
3638	nvan(:) = nvan_usda(njsc(:))
3639	avan(:) = avan_usda(njsc(:))
3640	mcr(:) = mcr_usda(njsc(:))
3641	mcs(:) = mcs_usda(njsc(:))
3642	ks(:) = ks_usda(njsc(:))
3643	!!$ mcfc(:) = mcf_usda(njsc(:))
3644	!!$ mcw(:) = mcw_usda(njsc(:))
3645
3646	!! Calculation of FC and WP based on above 5 parameters
3647	mvan(:) = un - (un / nvan(:))
3648	! Define matrix potential in mm for wilting point and field capacity (with sand vs clay-silt variation)
3649	psi_w(:) = 150000.
3650	DO ib=1, nbpt
3651	IF ( ks(ib) .GE. 560 ) THEN ! Sandy soils (560 is equivalent of 2.75 at log scale of Ks, mm/d)
3652	psi_fc(ib) = 1000.
3653	ELSE ! Finer soils
3654	psi_fc(ib) = 3300.
3655	ENDIF
3656	ENDDO
3657	mcfc(:) = mcr(:) + (( mcs(:) - mcr(:)) / (un + ( avan(:) * psi_fc(:)) nvan(:)) mvan(:))
3658	mcw(:) = mcr(:) + (( mcs(:) - mcr(:)) / (un + ( avan(:) * psi_w(:)) nvan(:)) mvan(:))
3659
3660	ELSE ! if spmipexp is not among texture or maps or unif
3661	WRITE(numout,*) "Unsupported spmipexp=",spmipexp
3662	WRITE(numout,*) "Choose between texture, maps, and unif"
3663	CALL ipslerr_p(3,'soilproc_soilt','Bad choice of spmipexp','Choose between texture, maps, and unif','')
3664	ENDIF
3665	ENDIF
3666	ENDIF ! SPMIPEXP
3667
3668	! Write diagnostics
3669	CALL xios_orchidee_send_field("atext",atext)
3670	CALL xios_orchidee_send_field("interp_diag_atext",atext)
3671	CALL xios_orchidee_send_field("interp_diag_soilclass",soilclass)
3672	CALL xios_orchidee_send_field("interp_diag_clayfraction",clayfraction)
3673
3674	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_soilt ended'
3675
3676	END SUBROUTINE slowproc_soilt
3677
3678	!! ================================================================================================================================
3679	!! SUBROUTINE : slowproc_slope
3680	!!
3681	!>\BRIEF Calculate mean slope coef in each model grid box from the slope map
3682	!!
3683	!! DESCRIPTION : (definitions, functional, design, flags):
3684	!!
3685	!! RECENT CHANGE(S): None
3686	!!
3687	!! MAIN OUTPUT VARIABLE(S): ::reinf_slope
3688	!!
3689	!! REFERENCE(S) : None
3690	!!
3691	!! FLOWCHART : None
3692	!! \n
3693	!_ ================================================================================================================================
3694
3695	SUBROUTINE slowproc_slope(nbpt, lalo, neighbours, resolution, contfrac, reinf_slope)
3696
3697	USE interpweight
3698
3699	IMPLICIT NONE
3700
3701	!
3702	!
3703	!
3704	! 0.1 INPUT
3705	!
3706	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3707	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
3708	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point
3709	! (1=North and then clockwise)
3710	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
3711	REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Fraction of continent in the grid
3712	!
3713	! 0.2 OUTPUT
3714	!
3715	REAL(r_std), INTENT(out) :: reinf_slope(nbpt) ! slope coef
3716	!
3717	! 0.3 LOCAL
3718	!
3719	!
3720	REAL(r_std) :: slope_noreinf ! Slope above which runoff is maximum
3721	CHARACTER(LEN=80) :: filename
3722	REAL(r_std) :: vmin, vmax !! min/max values to use for the
3723	!! renormalization
3724	REAL(r_std), DIMENSION(nbpt) :: aslope !! slope availability
3725
3726	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3727	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
3728	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3729	!! 'XYKindTime': Input values are kinds
3730	!! of something with a temporal
3731	!! evolution on the dx*dy matrix'
3732	LOGICAL :: nonegative !! whether negative values should be removed
3733	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3734	!! 'nomask': no-mask is applied
3735	!! 'mbelow': take values below maskvals(1)
3736	!! 'mabove': take values above maskvals(1)
3737	!! 'msumrange': take values within 2 ranges;
3738	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3739	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3740	!! (normalized by maskvals(3))
3741	!! 'var': mask values are taken from a
3742	!! variable inside the file (>0)
3743	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3744	!! `maskingtype')
3745	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3746
3747	!_ ================================================================================================================================
3748
3749	!
3750	!Config Key = SLOPE_NOREINF
3751	!Config Desc = Slope over which surface runoff does not reinfiltrate
3752	!Config If =
3753	!Config Def = 0.5
3754	!Config Help = The slope above which there is no reinfiltration
3755	!Config Units = [%]
3756	!
3757	slope_noreinf = 0.5 ! slope in percent
3758	!
3759	CALL getin_p('SLOPE_NOREINF',slope_noreinf)
3760	!
3761	!Config Key = TOPOGRAPHY_FILE
3762	!Config Desc = Name of file from which the topography map is to be read
3763	!Config If =
3764	!Config Def = cartepente2d_15min.nc
3765	!Config Help = The name of the file to be opened to read the orography
3766	!Config map is to be given here. Usualy SECHIBA runs with a 2'
3767	!Config map which is derived from the NGDC one.
3768	!Config Units = [FILE]
3769	!
3770	filename = 'cartepente2d_15min.nc'
3771	CALL getin_p('TOPOGRAPHY_FILE',filename)
3772
3773	IF (xios_interpolation) THEN
3774
3775	CALL xios_orchidee_recv_field('reinf_slope_interp',reinf_slope)
3776	CALL xios_orchidee_recv_field('frac_slope_interp',aslope)
3777
3778
3779	ELSE
3780
3781	variablename = 'pente'
3782	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_slope: Read and interpolate " &
3783	// TRIM(filename) // " for variable " // TRIM(variablename)
3784
3785	! For this case there are not types/categories. We have 'only' a continuos field
3786	! Assigning values to vmin, vmax
3787	vmin = 0.
3788	vmax = 9999.
3789
3790	!! Variables for interpweight
3791	! Type of calculation of cell fractions
3792	fractype = 'slopecalc'
3793	! Name of the longitude and latitude in the input file
3794	lonname = 'longitude'
3795	latname = 'latitude'
3796	! Should negative values be set to zero from input file?
3797	nonegative = .FALSE.
3798	! Type of mask to apply to the input data (see header for more details)
3799	maskingtype = 'mabove'
3800	! Values to use for the masking
3801	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3802	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
3803	namemaskvar = ''
3804
3805	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3806	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3807	maskvals, namemaskvar, -1, fractype, slope_default, slope_noreinf, &
3808	reinf_slope, aslope)
3809	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_slope after interpweight_2Dcont'
3810
3811	ENDIF
3812
3813	! Write diagnostics
3814	CALL xios_orchidee_send_field("aslope",aslope)
3815	CALL xios_orchidee_send_field("interp_diag_aslope",aslope)
3816
3817	CALL xios_orchidee_send_field("interp_diag_reinf_slope",reinf_slope)
3818
3819	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_slope ended'
3820
3821	END SUBROUTINE slowproc_slope
3822
3823
3824	!! ================================================================================================================================
3825	!! SUBROUTINE : slowproc_woodharvest
3826	!!
3827	!>\BRIEF
3828	!!
3829	!! DESCRIPTION :
3830	!!
3831	!! RECENT CHANGE(S): None
3832	!!
3833	!! MAIN OUTPUT VARIABLE(S): ::
3834	!!
3835	!! REFERENCE(S) : None
3836	!!
3837	!! FLOWCHART : None
3838	!! \n
3839	!_ ================================================================================================================================
3840
3841	SUBROUTINE slowproc_woodharvest(nbpt, lalo, neighbours, resolution, contfrac, woodharvest)
3842
3843	USE interpweight
3844
3845	IMPLICIT NONE
3846
3847	!
3848	!
3849	!
3850	! 0.1 INPUT
3851	!
3852	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
3853	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
3854	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
3855	!! (1=North and then clockwise)
3856	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
3857	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
3858	!
3859	! 0.2 OUTPUT
3860	!
3861	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: woodharvest !! Wood harvest
3862	!
3863	! 0.3 LOCAL
3864	!
3865	CHARACTER(LEN=80) :: filename
3866	REAL(r_std) :: vmin, vmax
3867	REAL(r_std), DIMENSION(nbpt) :: aoutvar !! availability of input data to
3868	!! interpolate output variable
3869	!! (on the nbpt space)
3870	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3871	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
3872	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3873	!! 'XYKindTime': Input values are kinds
3874	!! of something with a temporal
3875	!! evolution on the dx*dy matrix'
3876	LOGICAL :: nonegative !! whether negative values should be removed
3877	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3878	!! 'nomask': no-mask is applied
3879	!! 'mbelow': take values below maskvals(1)
3880	!! 'mabove': take values above maskvals(1)
3881	!! 'msumrange': take values within 2 ranges;
3882	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3883	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3884	!! (normalized by maskvals(3))
3885	!! 'var': mask values are taken from a
3886	!! variable inside the file (>0)
3887	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3888	!! `maskingtype')
3889	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3890	REAL(r_std), DIMENSION(1) :: variabletypevals !!
3891	! REAL(r_std), DIMENSION(nbp_mpi) :: woodharvest_mpi !! Wood harvest where all thredds OMP are gatherd
3892	!_ ================================================================================================================================
3893
3894
3895	!Config Key = WOODHARVEST_FILE
3896	!Config Desc = Name of file from which the wood harvest will be read
3897	!Config If = DO_WOOD_HARVEST
3898	!Config Def = woodharvest.nc
3899	!Config Help =
3900	!Config Units = [FILE]
3901	filename = 'woodharvest.nc'
3902	CALL getin_p('WOODHARVEST_FILE',filename)
3903	variablename = 'woodharvest'
3904
3905
3906	IF (xios_interpolation) THEN
3907	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readwoodharvest: Use XIOS to read and interpolate " &
3908	// TRIM(filename) // " for variable " // TRIM(variablename)
3909
3910	CALL xios_orchidee_recv_field('woodharvest_interp',woodharvest)
3911
3912	aoutvar = 1.0
3913	ELSE
3914
3915	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readwoodharvest: Read and interpolate " &
3916	// TRIM(filename) // " for variable " // TRIM(variablename)
3917
3918	! For this case there are not types/categories. We have 'only' a continuos field
3919	! Assigning values to vmin, vmax
3920	vmin = 0.
3921	vmax = 9999.
3922
3923	!! Variables for interpweight
3924	! Type of calculation of cell fractions
3925	fractype = 'default'
3926	! Name of the longitude and latitude in the input file
3927	lonname = 'longitude'
3928	latname = 'latitude'
3929	! Should negative values be set to zero from input file?
3930	nonegative = .TRUE.
3931	! Type of mask to apply to the input data (see header for more details)
3932	maskingtype = 'nomask'
3933	! Values to use for the masking
3934	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3935	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
3936	namemaskvar = ''
3937
3938	variabletypevals=-un
3939	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3940	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3941	maskvals, namemaskvar, -1, fractype, 0., 0., woodharvest, aoutvar)
3942	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_wodharvest after interpweight_2Dcont'
3943
3944	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_woodharvest ended'
3945	END IF
3946	END SUBROUTINE slowproc_woodharvest
3947
3948
3949	!! ================================================================================================================================
3950	!! SUBROUTINE : get_soilcorr_zobler
3951	!!
3952	!>\BRIEF The "get_soilcorr" routine defines the table of correspondence
3953	!! between the Zobler types and the three texture types known by SECHIBA and STOMATE :
3954	!! silt, sand and clay.
3955	!!
3956	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
3957	!! The data from this file is then interpolated to the grid of the model. \n
3958	!! The aim is to get fractions for sand loam and clay in each grid box.\n
3959	!! This information is used for soil hydrology and respiration.
3960	!!
3961	!!
3962	!! RECENT CHANGE(S): None
3963	!!
3964	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
3965	!!
3966	!! REFERENCE(S) :
3967	!! - Zobler L., 1986, A World Soil File for global climate modelling. NASA Technical memorandum 87802. NASA
3968	!! Goddard Institute for Space Studies, New York, U.S.A.
3969	!!
3970	!! FLOWCHART : None
3971	!! \n
3972	!_ ================================================================================================================================
3973
3974	SUBROUTINE get_soilcorr_zobler (nzobler,textfrac_table)
3975
3976	IMPLICIT NONE
3977
3978	!! 0. Variables and parameters declaration
3979
3980	INTEGER(i_std),PARAMETER :: nbtypes_zobler = 7 !! Number of Zobler types (unitless)
3981
3982	!! 0.1 Input variables
3983
3984	INTEGER(i_std),INTENT(in) :: nzobler !! Size of the array (unitless)
3985
3986	!! 0.2 Output variables
3987
3988	REAL(r_std),DIMENSION(nzobler,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
3989	!! and granulometric composition (0-1, unitless)
3990
3991	!! 0.4 Local variables
3992
3993	INTEGER(i_std) :: ib !! Indice (unitless)
3994
3995	!_ ================================================================================================================================
3996
3997	!-
3998	! 0. Check consistency
3999	!-
4000	IF (nzobler /= nbtypes_zobler) THEN
4001	CALL ipslerr_p(3,'get_soilcorr', 'nzobler /= nbtypes_zobler',&
4002	& 'We do not have the correct number of classes', &
4003	& ' in the code for the file.') ! Fatal error
4004	ENDIF
4005
4006	!-
4007	! 1. Textural fraction for : silt sand clay
4008	!-
4009	textfrac_table(1,:) = (/ 0.12, 0.82, 0.06 /)
4010	textfrac_table(2,:) = (/ 0.32, 0.58, 0.10 /)
4011	textfrac_table(3,:) = (/ 0.39, 0.43, 0.18 /)
4012	textfrac_table(4,:) = (/ 0.15, 0.58, 0.27 /)
4013	textfrac_table(5,:) = (/ 0.34, 0.32, 0.34 /)
4014	textfrac_table(6,:) = (/ 0.00, 1.00, 0.00 /)
4015	textfrac_table(7,:) = (/ 0.39, 0.43, 0.18 /)
4016
4017
4018	!-
4019	! 2. Check the mapping for the Zobler types which are going into the ORCHIDEE textures classes
4020	!-
4021	DO ib=1,nzobler ! Loop over # classes soil
4022
4023	IF (ABS(SUM(textfrac_table(ib,:))-1.0) > EPSILON(1.0)) THEN ! The sum of the textural fractions should not exceed 1 !
4024	WRITE(numout,*) &
4025	& 'Error in the correspondence table', &
4026	& ' sum is not equal to 1 in', ib
4027	WRITE(numout,*) textfrac_table(ib,:)
4028	CALL ipslerr_p(3,'get_soilcorr', 'SUM(textfrac_table(ib,:)) /= 1.0',&
4029	& '', 'Error in the correspondence table') ! Fatal error
4030	ENDIF
4031
4032	ENDDO ! Loop over # classes soil
4033
4034
4035	END SUBROUTINE get_soilcorr_zobler
4036
4037	!! ================================================================================================================================
4038	!! SUBROUTINE : get_soilcorr_usda
4039	!!
4040	!>\BRIEF The "get_soilcorr_usda" routine defines the table of correspondence
4041	!! between the 12 USDA textural classes and their granulometric composition,
4042	!! as % of silt, sand and clay. This is used to further defien clayfraction.
4043	!!
4044	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
4045	!! The data from this file is then interpolated to the grid of the model. \n
4046	!! The aim is to get fractions for sand loam and clay in each grid box.\n
4047	!! This information is used for soil hydrology and respiration.
4048	!! The default map in this case is derived from Reynolds et al 2000, \n
4049	!! at the 1/12deg resolution, with indices that are consistent with the \n
4050	!! textures tabulated below
4051	!!
4052	!! RECENT CHANGE(S): Created by A. Ducharne on July 02, 2014
4053	!!
4054	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
4055	!!
4056	!! REFERENCE(S) :
4057	!!
4058	!! FLOWCHART : None
4059	!! \n
4060	!_ ================================================================================================================================
4061
4062	SUBROUTINE get_soilcorr_usda (nusda,textfrac_table)
4063
4064	IMPLICIT NONE
4065
4066	!! 0. Variables and parameters declaration
4067
4068	!! 0.1 Input variables
4069
4070	INTEGER(i_std),INTENT(in) :: nusda !! Size of the array (unitless)
4071
4072	!! 0.2 Output variables
4073
4074	REAL(r_std),DIMENSION(nusda,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
4075	!! and granulometric composition (0-1, unitless)
4076
4077	!! 0.4 Local variables
4078
4079	INTEGER(i_std),PARAMETER :: nbtypes_usda = 13 !! Number of USDA texture classes (unitless)
4080	INTEGER(i_std) :: n !! Index (unitless)
4081
4082	!_ ================================================================================================================================
4083
4084	!-
4085	! 0. Check consistency
4086	!-
4087	IF (nusda /= nbtypes_usda) THEN
4088	CALL ipslerr_p(3,'get_soilcorr', 'nusda /= nbtypes_usda',&
4089	& 'We do not have the correct number of classes', &
4090	& ' in the code for the file.') ! Fatal error
4091	ENDIF
4092
4093	!! Parameters for soil type distribution :
4094	!! Sand, Loamy Sand, Sandy Loam, Silt Loam, Silt, Loam, Sandy Clay Loam, Silty Clay Loam, Clay Loam, Sandy Clay, Silty Clay, Clay
4095	! The order comes from constantes_soil.f90
4096	! The corresponding granulometric composition comes from Carsel & Parrish, 1988
4097
4098	!-
4099	! 1. Textural fractions for : sand, clay
4100	!-
4101	textfrac_table(1,2:3) = (/ 0.93, 0.03 /) ! Sand
4102	textfrac_table(2,2:3) = (/ 0.81, 0.06 /) ! Loamy Sand
4103	textfrac_table(3,2:3) = (/ 0.63, 0.11 /) ! Sandy Loam
4104	textfrac_table(4,2:3) = (/ 0.17, 0.19 /) ! Silt Loam
4105	textfrac_table(5,2:3) = (/ 0.06, 0.10 /) ! Silt
4106	textfrac_table(6,2:3) = (/ 0.40, 0.20 /) ! Loam
4107	textfrac_table(7,2:3) = (/ 0.54, 0.27 /) ! Sandy Clay Loam
4108	textfrac_table(8,2:3) = (/ 0.08, 0.33 /) ! Silty Clay Loam
4109	textfrac_table(9,2:3) = (/ 0.30, 0.33 /) ! Clay Loam
4110	textfrac_table(10,2:3) = (/ 0.48, 0.41 /) ! Sandy Clay
4111	textfrac_table(11,2:3) = (/ 0.06, 0.46 /) ! Silty Clay
4112	textfrac_table(12,2:3) = (/ 0.15, 0.55 /) ! Clay
4113	textfrac_table(13,2:3) = (/ 0.15, 0.55 /) ! Clay
4114
4115	! Fraction of silt
4116
4117	DO n=1,nusda
4118	textfrac_table(n,1) = 1. - textfrac_table(n,2) - textfrac_table(n,3)
4119	END DO
4120
4121	END SUBROUTINE get_soilcorr_usda
4122
4123	!! ================================================================================================================================
4124	!! FUNCTION : tempfunc
4125	!!
4126	!>\BRIEF ! This function interpolates value between ztempmin and ztempmax
4127	!! used for lai detection.
4128	!!
4129	!! DESCRIPTION : This subroutine calculates a scalar between 0 and 1 with the following equation :\n
4130	!! \latexonly
4131	!! \input{constantes_veg_tempfunc.tex}
4132	!! \endlatexonly
4133	!!
4134	!! RECENT CHANGE(S): None
4135	!!
4136	!! RETURN VALUE : tempfunc_result
4137	!!
4138	!! REFERENCE(S) : None
4139	!!
4140	!! FLOWCHART : None
4141	!! \n
4142	!_ ================================================================================================================================
4143
4144	FUNCTION tempfunc (temp_in) RESULT (tempfunc_result)
4145
4146
4147	!! 0. Variables and parameters declaration
4148
4149	REAL(r_std),PARAMETER :: ztempmin=273._r_std !! Temperature for laimin (K)
4150	REAL(r_std),PARAMETER :: ztempmax=293._r_std !! Temperature for laimax (K)
4151	REAL(r_std) :: zfacteur !! Interpolation factor (K^{-2})
4152
4153	!! 0.1 Input variables
4154
4155	REAL(r_std),INTENT(in) :: temp_in !! Temperature (K)
4156
4157	!! 0.2 Result
4158
4159	REAL(r_std) :: tempfunc_result !! (unitless)
4160
4161	!_ ================================================================================================================================
4162
4163	!! 1. Define a coefficient
4164	zfacteur = un/(ztempmax-ztempmin)**2
4165
4166	!! 2. Computes tempfunc
4167	IF (temp_in > ztempmax) THEN
4168	tempfunc_result = un
4169	ELSEIF (temp_in < ztempmin) THEN
4170	tempfunc_result = zero
4171	ELSE
4172	tempfunc_result = un-zfacteur(ztempmax-temp_in)*2
4173	ENDIF !(temp_in > ztempmax)
4174
4175
4176	END FUNCTION tempfunc
4177
4178
4179	!! ================================================================================================================================
4180	!! SUBROUTINE : slowproc_checkveget
4181	!!
4182	!>\BRIEF To verify the consistency of the various fractions defined within the grid box after having been
4183	!! been updated by STOMATE or the standard procedures.
4184	!!
4185	!! DESCRIPTION : (definitions, functional, design, flags):
4186	!!
4187	!! RECENT CHANGE(S): None
4188	!!
4189	!! MAIN OUTPUT VARIABLE(S): :: none
4190	!!
4191	!! REFERENCE(S) : None
4192	!!
4193	!! FLOWCHART : None
4194	!! \n
4195	!_ ================================================================================================================================
4196	!
4197	SUBROUTINE slowproc_checkveget(nbpt, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
4198
4199	! 0.1 INPUT
4200	!
4201	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
4202	REAL(r_std),DIMENSION (nbpt,nnobio), INTENT(in) :: frac_nobio ! Fraction of ice,lakes,cities, ... (unitless)
4203	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget_max ! Maximum fraction of vegetation type including none biological fraction (unitless)
4204	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget ! Vegetation fractions
4205	REAL(r_std),DIMENSION (nbpt), INTENT(in) :: tot_bare_soil ! Total evaporating bare soil fraction within the mesh
4206	REAL(r_std),DIMENSION (nbpt,nstm), INTENT(in) :: soiltile ! Fraction of soil tiles in the gridbox (unitless)
4207
4208	! 0.3 LOCAL
4209	!
4210	INTEGER(i_std) :: ji, jn, jv
4211	REAL(r_std) :: epsilocal !! A very small value
4212	REAL(r_std) :: totfrac
4213	CHARACTER(len=80) :: str1, str2
4214
4215	!_ ================================================================================================================================
4216
4217	!
4218	! There is some margin added as the computing errors might bring us above EPSILON(un)
4219	!
4220	epsilocal = EPSILON(un)*1000.
4221
4222	!! 1.0 Verify that none of the fractions are smaller than min_vegfrac, without beeing zero.
4223	!!
4224	DO ji=1,nbpt
4225	DO jn=1,nnobio
4226	IF ( frac_nobio(ji,jn) > epsilocal .AND. frac_nobio(ji,jn) < min_vegfrac ) THEN
4227	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
4228	WRITE(str2,'("The small value obtained is ", E14.4)') frac_nobio(ji,jn)
4229	CALL ipslerr_p (3,'slowproc_checkveget', &
4230	"frac_nobio is larger than zero but smaller than min_vegfrac.", str1, str2)
4231	ENDIF
4232	ENDDO
4233	END DO
4234
4235	IF (.NOT. ok_dgvm) THEN
4236	DO ji=1,nbpt
4237	DO jv=1,nvm
4238	IF ( veget_max(ji,jv) > epsilocal .AND. veget_max(ji,jv) < 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)') veget_max(ji,jv)
4241	CALL ipslerr_p (3,'slowproc_checkveget', &
4242	"veget_max is larger than zero but smaller than min_vegfrac.", str1, str2)
4243	ENDIF
4244	ENDDO
4245	ENDDO
4246	END IF
4247
4248	!! 2.0 verify that with all the fractions we cover the entire grid box
4249	!!
4250	DO ji=1,nbpt
4251	totfrac = zero
4252	DO jn=1,nnobio
4253	totfrac = totfrac + frac_nobio(ji,jn)
4254	ENDDO
4255	DO jv=1,nvm
4256	totfrac = totfrac + veget_max(ji,jv)
4257	ENDDO
4258	IF ( ABS(totfrac - un) > epsilocal) THEN
4259	WRITE(str1,'("This occurs on grid box", I8)') ji
4260	WRITE(str2,'("The sum over all fraction and error are ", E14.4, E14.4)') totfrac, ABS(totfrac - un)
4261	CALL ipslerr_p (3,'slowproc_checkveget', &
4262	"veget_max + frac_nobio is not equal to 1.", str1, str2)
4263	WRITE(,) "EPSILON =", epsilocal
4264	ENDIF
4265	ENDDO
4266
4267	!! 3.0 Verify that veget is smaller or equal to veget_max
4268	!!
4269	DO ji=1,nbpt
4270	DO jv=1,nvm
4271	IF ( jv == ibare_sechiba ) THEN
4272	IF ( ABS(veget(ji,jv) - veget_max(ji,jv)) > epsilocal ) THEN
4273	WRITE(str1,'("This occurs on grid box", I8)') ji
4274	WRITE(str2,'("The difference is ", E14.4)') veget(ji,jv) - veget_max(ji,jv)
4275	CALL ipslerr_p (3,'slowproc_checkveget', &
4276	"veget is not equal to veget_max on bare soil.", str1, str2)
4277	ENDIF
4278	ELSE
4279	IF ( veget(ji,jv) > veget_max(ji,jv) ) THEN
4280	WRITE(str1,'("This occurs on grid box", I8)') ji
4281	WRITE(str2,'("The values for veget and veget_max :", F8.4, F8.4)') veget(ji,jv), veget_max(ji,jv)
4282	CALL ipslerr_p (3,'slowproc_checkveget', &
4283	"veget is greater than veget_max.", str1, str2)
4284	ENDIF
4285	ENDIF
4286	ENDDO
4287	ENDDO
4288
4289	!! 4.0 Test tot_bare_soil in relation to the other variables
4290	!!
4291	DO ji=1,nbpt
4292	totfrac = zero
4293	DO jv=1,nvm
4294	totfrac = totfrac + (veget_max(ji,jv) - veget(ji,jv))
4295	ENDDO
4296	! add the bare soil fraction to totfrac
4297	totfrac = totfrac + veget(ji,ibare_sechiba)
4298	! do the test
4299	IF ( ABS(totfrac - tot_bare_soil(ji)) > epsilocal ) THEN
4300	WRITE(str1,'("This occurs on grid box", I8)') ji
4301	WRITE(str2,'("The values for tot_bare_soil, tot frac and error :", F8.4, F8.4, E14.4)') &
4302	& tot_bare_soil(ji), totfrac, ABS(totfrac - tot_bare_soil(ji))
4303	CALL ipslerr_p (3,'slowproc_checkveget', &
4304	"tot_bare_soil does not correspond to the total bare soil fraction.", str1, str2)
4305	ENDIF
4306	ENDDO
4307
4308	!! 5.0 Test that soiltile has the right sum
4309	!!
4310	DO ji=1,nbpt
4311	totfrac = SUM(soiltile(ji,:))
4312	IF ( ABS(totfrac - un) > epsilocal ) THEN
4313	WRITE(numout,*) "soiltile does not sum-up to one. This occurs on grid box", ji
4314	WRITE(numout,*) "The soiltile for ji are :", soiltile(ji,:)
4315	CALL ipslerr_p (2,'slowproc_checkveget', &
4316	"soiltile does not sum-up to one.", "", "")
4317	ENDIF
4318	ENDDO
4319
4320	END SUBROUTINE slowproc_checkveget
4321
4322
4323	!! ================================================================================================================================
4324	!! SUBROUTINE : slowproc_change_frac
4325	!!
4326	!>\BRIEF Update the vegetation fractions
4327	!!
4328	!! DESCRIPTION : Update the vegetation fractions. This subroutine is called in the same time step as lcchange in stomatelpj has
4329	!! has been done. This subroutine is called after the diagnostics have been written in sechiba_main.
4330	!!
4331	!! RECENT CHANGE(S): None
4332	!!
4333	!! MAIN OUTPUT VARIABLE(S): :: veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile
4334	!!
4335	!! REFERENCE(S) : None
4336	!!
4337	!! FLOWCHART : None
4338	!! \n
4339	!_ ================================================================================================================================
4340
4341	SUBROUTINE slowproc_change_frac(kjpindex, lai, &
4342	veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
4343	!
4344	! 0. Declarations
4345	!
4346	! 0.1 Input variables
4347	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
4348	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m^2 m^{-2})
4349
4350	! 0.2 Output variables
4351	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
4352	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget !! Fraction of vegetation type in the mesh (unitless)
4353	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
4354	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
4355	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
4356	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
4357	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
4358	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdfraclut !! Fraction of non woody vegetation in each landuse tile (0-1, unitless)
4359
4360	! 0.3 Local variables
4361	INTEGER(i_std) :: ji, jv !! Loop index
4362
4363
4364	!! Update vegetation fractions with the values coming from the vegetation file read in slowproc_readvegetmax.
4365	!! Partial update has been taken into account for the case with DGVM and AGRICULTURE in slowproc_readvegetmax.
4366	veget_max = veget_max_new
4367	frac_nobio = frac_nobio_new
4368
4369	!! Verification and correction on veget_max, calculation of veget and soiltile.
4370	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
4371
4372	!! Calculate tot_bare_soil needed in hydrol, diffuco and condveg (fraction of bare soil in the mesh)
4373	tot_bare_soil(:) = veget_max(:,1)
4374	DO jv = 2, nvm
4375	DO ji =1, kjpindex
4376	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
4377	ENDDO
4378	END DO
4379
4380	!! Do some basic tests on the surface fractions updated above
4381	CALL slowproc_checkveget(kjpindex, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
4382
4383	END SUBROUTINE slowproc_change_frac
4384
4385	END MODULE slowproc

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

Download in other formats: