Context Navigation

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

Last change on this file since 7442 was 7442, checked in by agnes.ducharne, 2 years ago
Corrrects units and improves comments for REINF_SLOPE and SLOPE_NOREINF.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 214.2 KB

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

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

Download in other formats: