Context Navigation

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

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

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

Download in other formats: