Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_sechiba/slowproc.f90 @ 7346

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

Update for ESM configuartion with carbon transport, P Cadule, JG :

4 optional variables are removed from the interface to LMDZ. These variables were not used in any official versions of LMDZ and therfor there are no problem with backward compatibility. The variables concerned are fields_in, fields_out in subroutine intersurf_initialize_gathered and fields_in_names, fields_out_names in subroutine intersurf_main_gathered.
Move parameter veget_update from slowproc into constantes_var

Property svn:keywords set to HeadURL Date Author Revision

File size: 214.1 KB

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

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

Download in other formats: