Context Navigation

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

Last change on this file since 5805 was 5605, checked in by josefine.ghattas, 6 years ago

Added new option VEGETMAP_RESET. This option should be used to to change vegetation map while keeping VEGET_UPDATE=0Y. This option will read the vegetation map even if it is stored in the restart file. It will also set to zero carbon related variables: prod10, prod100, flux10, flux100, convflux, cflux_prod10, cflux_prod100, convfluxpft(:,:)

Cadule, J. Ghattas

Property svn:keywords set to HeadURL Date Author Revision

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

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

Download in other formats: