Context Navigation

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

Last change on this file since 7502 was 7502, checked in by josefine.ghattas, 2 years ago
Add test if missmatch for land/sea mask in model and restart files. See ticket #806
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 214.9 KB

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

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

Download in other formats: