Context Navigation

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

Last change on this file since 6144 was 6144, checked in by josefine.ghattas, 5 years ago
As done in the trunk rev [6143]: Corrected bug on siltfraction when reading soils_param.nc file with option XIOS_INTERPOLATION=y and SOILTYPE_CLASSIF=zobler.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 182.6 KB

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

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

Download in other formats: