Context Navigation

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

Last change on this file since 6145 was 6145, checked in by josefine.ghattas, 5 years ago
Correction related to VEGET_UPDATE and small fractions as done in the trunk [6034]: Separated into new subroutine slowproc_veget_max_limit the code were veget_max and frac_nobio is set to 0 if smaller than min_vegfrac. This subroutine is now called in slowproc_veget (as before) and also called after the new map has been read in slowproc_readvegetmax. See ticket #456
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 184.9 KB

Line
1
2	! =================================================================================================================================
3	! MODULE : slowproc
4	!
5	! CONTACT : orchidee-help _at_ listes.ipsl.fr
6	!
7	! LICENCE : IPSL (2006)
8	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
9	!
10	!>\BRIEF Groups the subroutines that: (1) initialize all variables used in
11	!! slowproc_main, (2) prepare the restart file for the next simulation, (3) Update the
12	!! vegetation cover if needed, and (4) handle all slow processes if the carbon
13	!! cycle is activated (call STOMATE) or update the vegetation properties (LAI and
14	!! fractional cover) in the case of a run with only SECHIBA.
15	!!
16	!!\n DESCRIPTION: None
17	!!
18	!! RECENT CHANGE(S): Allowed reading of USDA map, Nov 2014, ADucharne
19	!!
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_max_limit
1659	!!
1660	!>\BRIEF Set small fractions of veget_max to zero and normalize to keep the sum equal 1
1661	!!
1662	!! DESCRIPTION : Set small fractions of veget_max to zero and normalize to keep the sum equal 1
1663	!!
1664	!! RECENT CHANGE(S): The subroutine was previously a part of slowproc_veget,
1665	!! but was separated to be called also from slowproc_readvegetmax in order
1666	!! to have limited/normalized vegetation fractions right after its reading
1667	!! from the file (added by V.Bastrikov, 15/06/2019)
1668	!!
1669	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, veget_max
1670	!!
1671	!! REFERENCE(S) : None
1672	!!
1673	!! FLOWCHART : None
1674	!! \n
1675	!_ ================================================================================================================================
1676
1677	SUBROUTINE slowproc_veget_max_limit (kjpindex, frac_nobio, veget_max)
1678	!
1679	! 0. Declarations
1680	!
1681	! 0.1 Input variables
1682	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1683
1684	! 0.2 Modified variables
1685	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
1686	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
1687
1688	! 0.4 Local scalar and varaiables
1689	INTEGER(i_std) :: ji, jv !! indices
1690	REAL(r_std) :: SUMveg !! Total vegetation summed across PFTs
1691
1692	!_ ================================================================================================================================
1693	IF (printlev_loc >= 3) WRITE(numout,*) 'Entering slowproc_veget_max_limit'
1694
1695	!! Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
1696	DO ji = 1, kjpindex
1697	IF ( SUM(frac_nobio(ji,:)) .LT. min_vegfrac ) THEN
1698	frac_nobio(ji,:) = zero
1699	ENDIF
1700
1701	IF (.NOT. ok_dgvm) THEN
1702	DO jv = 1, nvm
1703	IF ( veget_max(ji,jv) .LT. min_vegfrac ) THEN
1704	veget_max(ji,jv) = zero
1705	ENDIF
1706	ENDDO
1707	END IF
1708
1709	!! Normalize to keep the sum equal 1.
1710	SUMveg = SUM(frac_nobio(ji,:))+SUM(veget_max(ji,:))
1711	frac_nobio(ji,:) = frac_nobio(ji,:)/SUMveg
1712	veget_max(ji,:) = veget_max(ji,:)/SUMveg
1713	ENDDO
1714
1715	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_veget_max_limit ended'
1716
1717	END SUBROUTINE slowproc_veget_max_limit
1718
1719
1720	!! ================================================================================================================================
1721	!! SUBROUTINE : slowproc_veget
1722	!!
1723	!>\BRIEF Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
1724	!!
1725	!! DESCRIPTION : Set small fractions to zero and normalize to keep the sum equal 1. Calucate veget and soiltile.
1726	!! (1) Set veget_max and frac_nobio for fraction smaller than min_vegfrac.
1727	!! (2) Calculate veget
1728	!! (3) Calculate totfrac_nobio
1729	!! (4) Calculate soiltile
1730	!! (5) Calculate fraclut
1731	!!
1732	!! RECENT CHANGE(S): None
1733	!!
1734	!! MAIN OUTPUT VARIABLE(S): :: frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut
1735	!!
1736	!! REFERENCE(S) : None
1737	!!
1738	!! FLOWCHART : None
1739	!! \n
1740	!_ ================================================================================================================================
1741
1742	SUBROUTINE slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
1743	!
1744	! 0. Declarations
1745	!
1746	! 0.1 Input variables
1747	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1748	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(in) :: lai !! PFT leaf area index (m^{2} m^{-2})
1749
1750	! 0.2 Modified variables
1751	REAL(r_std), DIMENSION(kjpindex,nnobio), INTENT(inout) :: frac_nobio !! Fraction of the mesh which is covered by ice, lakes, ...
1752	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(inout) :: veget_max !! Maximum fraction of vegetation type including none biological fraction (unitless)
1753
1754	! 0.3 Output variables
1755	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: veget !! Fraction of pixel covered by PFT. Fraction accounts for none-biological land covers (unitless)
1756	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: totfrac_nobio
1757	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
1758	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
1759	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdFraclut !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
1760
1761	! 0.4 Local scalar and varaiables
1762	INTEGER(i_std) :: ji, jv, jst !! indices
1763
1764	!_ ================================================================================================================================
1765	IF (printlev_loc > 8) WRITE(numout,*) 'Entering slowproc_veget'
1766
1767	!! 1. Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
1768	!! Normalize to have the sum equal 1.
1769	CALL slowproc_veget_max_limit(kjpindex, frac_nobio, veget_max)
1770
1771	!! 2. Calculate veget
1772	!! If lai of a vegetation type (jv > 1) is small, increase soil part
1773	!! stomate-like calculation
1774	DO ji = 1, kjpindex
1775	veget(ji,1)=veget_max(ji,1)
1776	DO jv = 2, nvm
1777	veget(ji,jv) = veget_max(ji,jv) * ( un - exp( - lai(ji,jv) * ext_coeff_vegetfrac(jv) ) )
1778	ENDDO
1779	ENDDO
1780
1781
1782	!! 3. Calculate totfrac_nobio
1783	totfrac_nobio(:) = zero
1784	DO jv = 1, nnobio
1785	totfrac_nobio(:) = totfrac_nobio(:) + frac_nobio(:,jv)
1786	ENDDO
1787
1788
1789	!! 4. Calculate soiltiles
1790	!! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
1791	!! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
1792	!! in the other modules to perform separated energy balances
1793	! The sum of all soiltiles makes one, and corresponds to the bio fraction
1794	! of the grid cell (called vegtot in hydrol)
1795	soiltile(:,:) = zero
1796	DO jv = 1, nvm
1797	jst = pref_soil_veg(jv)
1798	DO ji = 1, kjpindex
1799	soiltile(ji,jst) = soiltile(ji,jst) + veget_max(ji,jv)
1800	ENDDO
1801	ENDDO
1802	DO ji = 1, kjpindex
1803	IF (totfrac_nobio(ji) .LT. (1-min_sechiba)) THEN
1804	soiltile(ji,:)=soiltile(ji,:)/(1.-totfrac_nobio(ji))
1805	ENDIF
1806	ENDDO
1807
1808	!! 5. Calculate fraction of landuse tiles to be used only for diagnostic variables
1809	fraclut(:,:)=0
1810	nwdFraclut(:,id_psl)=0
1811	nwdFraclut(:,id_crp)=1.
1812	nwdFraclut(:,id_urb)=xios_default_val
1813	nwdFraclut(:,id_pst)=xios_default_val
1814	DO jv=1,nvm
1815	IF (natural(jv)) THEN
1816	fraclut(:,id_psl) = fraclut(:,id_psl) + veget_max(:,jv)
1817	IF(.NOT. is_tree(jv)) THEN
1818	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl) + veget_max(:,jv)
1819	ENDIF
1820	ELSE
1821	fraclut(:,id_crp) = fraclut(:,id_crp) + veget_max(:,jv)
1822	ENDIF
1823	END DO
1824
1825	WHERE (fraclut(:,id_psl) > min_sechiba)
1826	nwdFraclut(:,id_psl) = nwdFraclut(:,id_psl)/fraclut(:,id_psl)
1827	ELSEWHERE
1828	nwdFraclut(:,id_psl) = xios_default_val
1829	END WHERE
1830
1831	END SUBROUTINE slowproc_veget
1832
1833
1834	!! ================================================================================================================================
1835	!! SUBROUTINE : slowproc_lai
1836	!!
1837	!>\BRIEF Do the interpolation of lai for the PFTs in case the laimap is not read
1838	!!
1839	!! DESCRIPTION : (definitions, functional, design, flags):
1840	!! (1) Interplation by using the mean value of laimin and laimax for the PFTs
1841	!! (2) Interpolation between laimax and laimin values by using the temporal
1842	!! variations
1843	!! (3) If problem occurs during the interpolation, the routine stops
1844	!!
1845	!! RECENT CHANGE(S): None
1846	!!
1847	!! MAIN OUTPUT VARIABLE(S): ::lai
1848	!!
1849	!! REFERENCE(S) : None
1850	!!
1851	!! FLOWCHART : None
1852	!! \n
1853	!_ ================================================================================================================================
1854
1855	SUBROUTINE slowproc_lai (kjpindex,lcanop,stempdiag,lalo,resolution,lai,laimap)
1856	!
1857	! 0. Declarations
1858	!
1859	!! 0.1 Input variables
1860	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
1861	INTEGER(i_std), INTENT(in) :: lcanop !! soil level used for LAI
1862	REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag !! Soil temperature (K) ???
1863	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
1864	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! Size in x an y of the grid (m) - surface area of the gridbox
1865	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: laimap !! map of lai read
1866
1867	!! 0.2 Output
1868	REAL(r_std), DIMENSION(kjpindex,nvm), INTENT(out) :: lai !! PFT leaf area index (m^{2} m^{-2})LAI
1869
1870	!! 0.4 Local
1871	INTEGER(i_std) :: ji,jv !! Local indices
1872	!_ ================================================================================================================================
1873
1874	!
1875	IF ( .NOT. read_lai ) THEN
1876
1877	lai(: ,1) = zero
1878	! On boucle sur 2,nvm au lieu de 1,nvm
1879	DO jv = 2,nvm
1880	SELECT CASE (type_of_lai(jv))
1881
1882	CASE ("mean ")
1883	!
1884	! 1. do the interpolation between laimax and laimin
1885	!
1886	lai(:,jv) = undemi * (llaimax(jv) + llaimin(jv))
1887	!
1888	CASE ("inter")
1889	!
1890	! 2. do the interpolation between laimax and laimin
1891	!
1892	DO ji = 1,kjpindex
1893	lai(ji,jv) = llaimin(jv) + tempfunc(stempdiag(ji,lcanop)) * (llaimax(jv) - llaimin(jv))
1894	ENDDO
1895	!
1896	CASE default
1897	!
1898	! 3. Problem
1899	!
1900	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
1901	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
1902	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=false','')
1903	END SELECT
1904
1905	ENDDO
1906	!
1907	ELSE
1908	lai(: ,1) = zero
1909	! On boucle sur 2,nvm au lieu de 1,nvm
1910	DO jv = 2,nvm
1911
1912	SELECT CASE (type_of_lai(jv))
1913
1914	CASE ("mean ")
1915	!
1916	! 1. force MAXVAL of laimap on lai on this PFT
1917	!
1918	DO ji = 1,kjpindex
1919	lai(ji,jv) = MAXVAL(laimap(ji,jv,:))
1920	ENDDO
1921	!
1922	CASE ("inter")
1923	!
1924	! 2. do the interpolation between laimax and laimin
1925	!
1926	!
1927	! If January
1928	!
1929	IF (month_end .EQ. 1 ) THEN
1930	IF (day_end .LE. 15) THEN
1931	lai(:,jv) = laimap(:,jv,12)(1-(day_end+15)/30.) + laimap(:,jv,1)((day_end+15)/30.)
1932	ELSE
1933	lai(:,jv) = laimap(:,jv,1)(1-(day_end-15)/30.) + laimap(:,jv,2)((day_end-15)/30.)
1934	ENDIF
1935	!
1936	! If December
1937	!
1938	ELSE IF (month_end .EQ. 12) THEN
1939	IF (day_end .LE. 15) THEN
1940	lai(:,jv) = laimap(:,jv,11)(1-(day_end+15)/30.) + laimap(:,jv,12)((day_end+15)/30.)
1941	ELSE
1942	lai(:,jv) = laimap(:,jv,12)(1-(day_end-15)/30.) + laimap(:,jv,1)((day_end-15)/30.)
1943	ENDIF
1944	!
1945	! ELSE
1946	!
1947	ELSE
1948	IF (day_end .LE. 15) THEN
1949	lai(:,jv) = laimap(:,jv,month_end-1)(1-(day_end+15)/30.) + laimap(:,jv,month_end)((day_end+15)/30.)
1950	ELSE
1951	lai(:,jv) = laimap(:,jv,month_end)(1-(day_end-15)/30.) + laimap(:,jv,month_end+1)((day_end-15)/30.)
1952	ENDIF
1953	ENDIF
1954	!
1955	CASE default
1956	!
1957	! 3. Problem
1958	!
1959	WRITE (numout,*) 'This kind of lai choice is not possible. '// &
1960	' We stop with type_of_lai ',jv,' = ', type_of_lai(jv)
1961	CALL ipslerr_p(3,'slowproc_lai','Bad value for type_of_lai','read_lai=true','')
1962	END SELECT
1963
1964	ENDDO
1965	ENDIF
1966
1967	END SUBROUTINE slowproc_lai
1968
1969	!! ================================================================================================================================
1970	!! SUBROUTINE : slowproc_interlai
1971	!!
1972	!>\BRIEF Interpolate the LAI map to the grid of the model
1973	!!
1974	!! DESCRIPTION : (definitions, functional, design, flags):
1975	!!
1976	!! RECENT CHANGE(S): None
1977	!!
1978	!! MAIN OUTPUT VARIABLE(S): ::laimap
1979	!!
1980	!! REFERENCE(S) : None
1981	!!
1982	!! FLOWCHART : None
1983	!! \n
1984	!_ ================================================================================================================================
1985
1986	SUBROUTINE slowproc_interlai(nbpt, lalo, resolution, neighbours, contfrac, laimap)
1987
1988	USE interpweight
1989
1990	IMPLICIT NONE
1991
1992	!
1993	!
1994	!
1995	! 0.1 INPUT
1996	!
1997	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
1998	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
1999	!! (beware of the order = 1 : latitude, 2 : longitude)
2000	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
2001	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2002	!! (1=North and then clockwise)
2003	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
2004	!
2005	! 0.2 OUTPUT
2006	!
2007	REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) !! lai read variable and re-dimensioned
2008	!
2009	! 0.3 LOCAL
2010	!
2011	CHARACTER(LEN=80) :: filename !! name of the LAI map read
2012	INTEGER(i_std) :: ib, ip, jp, it, jv
2013	REAL(r_std) :: lmax, lmin, ldelta
2014	LOGICAL :: renormelize_lai ! flag to force LAI renormelization
2015	INTEGER :: ier
2016
2017	REAL(r_std), DIMENSION(nbpt) :: alaimap !! availability of the lai interpolation
2018	INTEGER, DIMENSION(4) :: invardims
2019	REAL(r_std), DIMENSION(nbpt,nvm,12) :: lairefrac !! lai fractions re-dimensioned
2020	REAL(r_std), DIMENSION(nbpt,nvm,12) :: fraclaiinterp !! lai fractions re-dimensioned
2021	REAL(r_std), DIMENSION(:), ALLOCATABLE :: vmin, vmax !! min/max values to use for the
2022	!! renormalization
2023	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2024	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2025	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
2026	!! (variabletypevals(1) = -un, not used)
2027	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2028	!! 'XYKindTime': Input values are kinds
2029	!! of something with a temporal
2030	!! evolution on the dx*dy matrix'
2031	LOGICAL :: nonegative !! whether negative values should be removed
2032	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2033	!! 'nomask': no-mask is applied
2034	!! 'mbelow': take values below maskvals(1)
2035	!! 'mabove': take values above maskvals(1)
2036	!! 'msumrange': take values within 2 ranges;
2037	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2038	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2039	!! (normalized by maskvals(3))
2040	!! 'var': mask values are taken from a
2041	!! variable inside the file (>0)
2042	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2043	!! `maskingtype')
2044	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2045	!_ ================================================================================================================================
2046
2047	!
2048	!Config Key = LAI_FILE
2049	!Config Desc = Name of file from which the vegetation map is to be read
2050	!Config If = LAI_MAP
2051	!Config Def = lai2D.nc
2052	!Config Help = The name of the file to be opened to read the LAI
2053	!Config map is to be given here. Usualy SECHIBA runs with a 5kmx5km
2054	!Config map which is derived from a Nicolas VIOVY one.
2055	!Config Units = [FILE]
2056	!
2057	filename = 'lai2D.nc'
2058	CALL getin_p('LAI_FILE',filename)
2059	variablename = 'LAI'
2060
2061	IF (xios_interpolation) THEN
2062	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_interlai: Use XIOS to read and interpolate " &
2063	// TRIM(filename) //" for variable " //TRIM(variablename)
2064
2065	CALL xios_orchidee_recv_field('lai_interp',lairefrac)
2066	CALL xios_orchidee_recv_field('frac_lai_interp',fraclaiinterp)
2067	alaimap(:) = fraclaiinterp(:,1,1)
2068	ELSE
2069
2070	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_interlai: Start interpolate " &
2071	// TRIM(filename) //" for variable " //TRIM(variablename)
2072
2073	! invardims: shape of variable in input file to interpolate
2074	invardims = interpweight_get_var4dims_file(filename, variablename)
2075	! Check coherence of dimensions read from the file
2076	IF (invardims(4) /= 12) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of time dimension in input file for lai','','')
2077	IF (invardims(3) /= nvm) CALL ipslerr_p(3,'slowproc_interlai','Wrong dimension of PFT dimension in input file for lai','','')
2078
2079	ALLOCATE(vmin(nvm),stat=ier)
2080	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmin','','')
2081
2082	ALLOCATE(vmax(nvm), STAT=ier)
2083	IF (ier /= 0) CALL ipslerr_p(3,'slowproc_interlai','Problem in allocation of variable vmax','','')
2084
2085
2086	! Assigning values to vmin, vmax
2087	vmin = un
2088	vmax = nvm*un
2089
2090	variabletypevals = -un
2091
2092	!! Variables for interpweight
2093	! Type of calculation of cell fractions
2094	fractype = 'default'
2095	! Name of the longitude and latitude in the input file
2096	lonname = 'longitude'
2097	latname = 'latitude'
2098	! Should negative values be set to zero from input file?
2099	nonegative = .TRUE.
2100	! Type of mask to apply to the input data (see header for more details)
2101	maskingtype = 'mbelow'
2102	! Values to use for the masking
2103	maskvals = (/ 20., undef_sechiba, undef_sechiba /)
2104	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2105	namemaskvar = ''
2106
2107	CALL interpweight_4D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
2108	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2109	maskvals, namemaskvar, nvm, invardims(4), -1, fractype, &
2110	-1., -1., lairefrac, alaimap)
2111
2112	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_interlai after interpweight_4D'
2113
2114	ENDIF
2115
2116
2117
2118	!
2119	!
2120	!Config Key = RENORM_LAI
2121	!Config Desc = flag to force LAI renormelization
2122	!Config If = LAI_MAP
2123	!Config Def = n
2124	!Config Help = If true, the laimap will be renormalize between llaimin and llaimax parameters.
2125	!Config Units = [FLAG]
2126	!
2127	renormelize_lai = .FALSE.
2128	CALL getin_p('RENORM_LAI',renormelize_lai)
2129
2130	!
2131	laimap(:,:,:) = zero
2132	!
2133	IF (printlev_loc >= 5) THEN
2134	WRITE(numout,*)' slowproc_interlai before starting loop nbpt:', nbpt
2135	END IF
2136
2137	! Assigning the right values and giving a value where information was not found
2138	DO ib=1,nbpt
2139	IF (alaimap(ib) < min_sechiba) THEN
2140	DO jv=1,nvm
2141	laimap(ib,jv,:) = (llaimax(jv)+llaimin(jv))/deux
2142	ENDDO
2143	ELSE
2144	DO jv=1, nvm
2145	DO it=1, 12
2146	laimap(ib,jv,it) = lairefrac(ib,jv,it)
2147	ENDDO
2148	ENDDO
2149	END IF
2150	ENDDO
2151	!
2152	! Normelize the read LAI by the values SECHIBA is used to
2153	!
2154	IF ( renormelize_lai ) THEN
2155	DO ib=1,nbpt
2156	DO jv=1, nvm
2157	lmax = MAXVAL(laimap(ib,jv,:))
2158	lmin = MINVAL(laimap(ib,jv,:))
2159	ldelta = lmax-lmin
2160	IF ( ldelta < min_sechiba) THEN
2161	! LAI constante ... keep it constant
2162	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)+(llaimax(jv)+llaimin(jv))/deux
2163	ELSE
2164	laimap(ib,jv,:) = (laimap(ib,jv,:)-lmin)/(lmax-lmin)*(llaimax(jv)-llaimin(jv))+llaimin(jv)
2165	ENDIF
2166	ENDDO
2167	ENDDO
2168	ENDIF
2169
2170	! Write diagnostics
2171	CALL xios_orchidee_send_field("alaimap",alaimap)
2172	CALL xios_orchidee_send_field("interp_diag_lai",laimap)
2173
2174	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_interlai ended'
2175
2176	END SUBROUTINE slowproc_interlai
2177
2178	!! ================================================================================================================================
2179	!! SUBROUTINE : slowproc_readvegetmax
2180	!!
2181	!>\BRIEF Read and interpolate a vegetation map (by pft)
2182	!!
2183	!! DESCRIPTION : (definitions, functional, design, flags):
2184	!!
2185	!! RECENT CHANGE(S): The subroutine was previously called slowproc_update.
2186	!!
2187	!! MAIN OUTPUT VARIABLE(S):
2188	!!
2189	!! REFERENCE(S) : None
2190	!!
2191	!! FLOWCHART : None
2192	!! \n
2193	!_ ================================================================================================================================
2194
2195	SUBROUTINE slowproc_readvegetmax(nbpt, lalo, neighbours, resolution, contfrac, veget_last, &
2196	veget_next, frac_nobio_next, veget_year, init)
2197
2198	USE interpweight
2199	IMPLICIT NONE
2200
2201	!
2202	!
2203	!
2204	! 0.1 INPUT
2205	!
2206	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs
2207	!! to be interpolated
2208	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
2209	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
2210	!! (1=North and then clockwise)
2211	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
2212	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
2213	!
2214	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(in) :: veget_last !! old max vegetfrac
2215	INTEGER(i_std), INTENT(in) :: veget_year !! first year for landuse (0 == NO TIME AXIS)
2216	LOGICAL, INTENT(in) :: init !! initialisation : in case of dgvm, it forces update of all PFTs
2217	!
2218	! 0.2 OUTPUT
2219	!
2220	REAL(r_std), DIMENSION(nbpt,nvm), INTENT(out) :: veget_next !! new max vegetfrac
2221	REAL(r_std), DIMENSION(nbpt,nnobio), INTENT(out) :: frac_nobio_next !! new fraction of the mesh which is
2222	!! covered by ice, lakes, ...
2223
2224	!
2225	! 0.3 LOCAL
2226	!
2227	!
2228	CHARACTER(LEN=80) :: filename
2229	INTEGER(i_std) :: ib, inobio, jv
2230	REAL(r_std) :: sumf, err, norm
2231	!
2232	! for DGVM case :
2233	REAL(r_std) :: sum_veg ! sum of vegets
2234	REAL(r_std) :: sum_nobio ! sum of nobios
2235	REAL(r_std) :: sumvAnthro_old, sumvAnthro ! last an new sum of antrhopic vegets
2236	REAL(r_std) :: rapport ! (S-B) / (S-A)
2237	LOGICAL :: partial_update ! if TRUE, partialy update PFT (only anthropic ones)
2238	! e.g. in case of DGVM and not init (optional parameter)
2239	REAL(r_std), DIMENSION(nbpt,nvm) :: vegetrefrac !! veget fractions re-dimensioned
2240	REAL(r_std), DIMENSION(nbpt) :: aveget !! Availability of the soilcol interpolation
2241	REAL(r_std), DIMENSION(nbpt,nvm) :: aveget_nvm !! Availability of the soilcol interpolation
2242	REAL(r_std), DIMENSION(nvm) :: vmin, vmax !! min/max values to use for the renormalization
2243	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2244	CHARACTER(LEN=80) :: lonname, latname !! lon, lat names in input file
2245	REAL(r_std), DIMENSION(nvm) :: variabletypevals !! Values for all the types of the variable
2246	!! (variabletypevals(1) = -un, not used)
2247	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2248	!! 'XYKindTime': Input values are kinds
2249	!! of something with a temporal
2250	!! evolution on the dx*dy matrix'
2251	LOGICAL :: nonegative !! whether negative values should be removed
2252	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2253	!! 'nomask': no-mask is applied
2254	!! 'mbelow': take values below maskvals(1)
2255	!! 'mabove': take values above maskvals(1)
2256	!! 'msumrange': take values within 2 ranges;
2257	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2258	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2259	!! (normalized by maskvals(3))
2260	!! 'var': mask values are taken from a
2261	!! variable inside the file (>0)
2262	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2263	!! `maskingtype')
2264	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2265	CHARACTER(LEN=250) :: msg
2266
2267	!_ ================================================================================================================================
2268
2269	IF (printlev_loc >= 5) PRINT *,' In slowproc_readvegetmax'
2270
2271	!
2272	!Config Key = VEGETATION_FILE
2273	!Config Desc = Name of file from which the vegetation map is to be read
2274	!Config If =
2275	!Config Def = PFTmap.nc
2276	!Config Help = The name of the file to be opened to read a vegetation
2277	!Config map (in pft) is to be given here.
2278	!Config Units = [FILE]
2279	!
2280	filename = 'PFTmap.nc'
2281	CALL getin_p('VEGETATION_FILE',filename)
2282	variablename = 'maxvegetfrac'
2283
2284
2285	IF (xios_interpolation) THEN
2286	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readvegetmax: Use XIOS to read and interpolate " &
2287	// TRIM(filename) // " for variable " // TRIM(variablename)
2288
2289	CALL xios_orchidee_recv_field('frac_veget',vegetrefrac)
2290	CALL xios_orchidee_recv_field('frac_veget_frac',aveget_nvm)
2291	aveget(:)=aveget_nvm(:,1)
2292
2293	DO ib = 1, nbpt
2294	IF (aveget(ib) > min_sechiba) THEN
2295	vegetrefrac(ib,:) = vegetrefrac(ib,:)/aveget(ib) ! intersected area normalization
2296	vegetrefrac(ib,:) = vegetrefrac(ib,:)/SUM(vegetrefrac(ib,:))
2297	ENDIF
2298	ENDDO
2299
2300	ELSE
2301
2302	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_readvegetmax: Start interpolate " &
2303	// TRIM(filename) // " for variable " // TRIM(variablename)
2304
2305	! Assigning values to vmin, vmax
2306	vmin = 1
2307	vmax = nvm*1._r_std
2308
2309	variabletypevals = -un
2310
2311	!! Variables for interpweight
2312	! Type of calculation of cell fractions
2313	fractype = 'default'
2314	! Name of the longitude and latitude in the input file
2315	lonname = 'lon'
2316	latname = 'lat'
2317	! Should negative values be set to zero from input file?
2318	nonegative = .FALSE.
2319	! Type of mask to apply to the input data (see header for more details)
2320	maskingtype = 'msumrange'
2321	! Values to use for the masking
2322	maskvals = (/ 1.-1.e-7, 0., 2. /)
2323	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
2324	namemaskvar = ''
2325
2326	CALL interpweight_3D(nbpt, nvm, variabletypevals, lalo, resolution, neighbours, &
2327	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2328	maskvals, namemaskvar, nvm, 0, veget_year, fractype, &
2329	-1., -1., vegetrefrac, aveget)
2330	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after interpeeight_3D'
2331	ENDIF
2332	!
2333	! Compute the logical for partial (only anthropic) PTFs update
2334	IF (ok_dgvm .AND. .NOT. init) THEN
2335	partial_update= .TRUE.
2336	ELSE
2337	partial_update=.FALSE.
2338	END IF
2339
2340	IF (printlev_loc >= 5) THEN
2341	WRITE(numout,*)' slowproc_readvegetmax before updating loop nbpt:', nbpt
2342	END IF
2343
2344	IF ( .NOT. partial_update ) THEN
2345	veget_next(:,:)=zero
2346
2347	IF (printlev_loc >=3 .AND. ANY(aveget < min_sechiba)) THEN
2348	WRITE(numout,*) 'Some grid cells on the model grid did not have any points on the source grid.'
2349	IF (init) THEN
2350	WRITE(numout,*) 'Initialization with full fraction of bare soil are done for the below grid cells.'
2351	ELSE
2352	WRITE(numout,*) 'Old values are kept for the below grid cells.'
2353	ENDIF
2354	WRITE(numout,*) 'List of grid cells (ib, lat, lon):'
2355	END IF
2356
2357	DO ib = 1, nbpt
2358	! vegetrefrac is already normalized to sum equal one for each grid cell
2359	veget_next(ib,:) = vegetrefrac(ib,:)
2360
2361	IF (aveget(ib) < min_sechiba) THEN
2362	IF (printlev_loc >=3) WRITE(numout,*) ib,lalo(ib,1),lalo(ib,2)
2363	IF (init) THEN
2364	veget_next(ib,1) = un
2365	veget_next(ib,2:nvm) = zero
2366	ELSE
2367	veget_next(ib,:) = veget_last(ib,:)
2368	ENDIF
2369	ENDIF
2370	ENDDO
2371	ELSE
2372	! Partial update
2373	DO ib = 1, nbpt
2374	IF (aveget(ib) > min_sechiba) THEN
2375	! For the case with properly interpolated grid cells (aveget>0)
2376
2377	! last veget for this point
2378	sum_veg=SUM(veget_last(ib,:))
2379	!
2380	! If the DGVM is activated, only anthropic PFTs are utpdated, the others are copied from previous time-step
2381	veget_next(ib,:) = veget_last(ib,:)
2382
2383	DO jv = 2, nvm
2384	IF ( .NOT. natural(jv) ) THEN
2385	veget_next(ib,jv) = vegetrefrac(ib,jv)
2386	ENDIF
2387	ENDDO
2388
2389	sumvAnthro_old = zero
2390	sumvAnthro = zero
2391	DO jv = 2, nvm
2392	IF ( .NOT. natural(jv) ) THEN
2393	sumvAnthro = sumvAnthro + veget_next(ib,jv)
2394	sumvAnthro_old = sumvAnthro_old + veget_last(ib,jv)
2395	ENDIF
2396	ENDDO
2397
2398	IF ( sumvAnthro_old < sumvAnthro ) THEN
2399	! Increase of non natural vegetations (increase of agriculture)
2400	! The proportion of natural PFT's must be preserved
2401	! ie the sum of vegets is preserved
2402	! and natural PFT / (sum of veget - sum of antropic veget)
2403	! is preserved.
2404	rapport = ( sum_veg - sumvAnthro ) / ( sum_veg - sumvAnthro_old )
2405	DO jv = 1, nvm
2406	IF ( natural(jv) ) THEN
2407	veget_next(ib,jv) = veget_last(ib,jv) * rapport
2408	ENDIF
2409	ENDDO
2410	ELSE
2411	! Increase of natural vegetations (decrease of agriculture)
2412	! The decrease of agriculture is replaced by bare soil. The DGVM will
2413	! re-introduce natural PFT's.
2414	DO jv = 1, nvm
2415	IF ( natural(jv) ) THEN
2416	veget_next(ib,jv) = veget_last(ib,jv)
2417	ENDIF
2418	ENDDO
2419	veget_next(ib,1) = veget_next(ib,1) + sumvAnthro_old - sumvAnthro
2420	ENDIF
2421
2422	! test
2423	IF ( ABS( SUM(veget_next(ib,:)) - sum_veg ) > 10*EPSILON(un) ) THEN
2424	WRITE(numout,*) 'slowproc_readvegetmax _______'
2425	msg = " No conservation of sum of veget for point "
2426	WRITE(numout,*) TRIM(msg), ib, ",(", lalo(ib,1),",", lalo(ib,2), ")"
2427	WRITE(numout,*) " last sum of veget ", sum_veg, " new sum of veget ", &
2428	SUM(veget_next(ib,:)), " error : ", SUM(veget_next(ib,:))-sum_veg
2429	WRITE(numout,*) " Anthropic modifications : last ",sumvAnthro_old," new ",sumvAnthro
2430	CALL ipslerr_p (3,'slowproc_readvegetmax', &
2431	& 'No conservation of sum of veget_next', &
2432	& "The sum of veget_next is different after reading Land Use map.", &
2433	& '(verify the dgvm case model.)')
2434	ENDIF
2435	ELSE
2436	! For the case when there was a propblem with the interpolation, aveget < min_sechiba
2437	WRITE(numout,*) 'slowproc_readvegetmax _______'
2438	WRITE(numout,*) " No land point in the map for point ", ib, ",(", lalo(ib,1), ",", &
2439	lalo(ib,2),")"
2440	CALL ipslerr_p (2,'slowproc_readvegetmax', &
2441	& 'Problem with vegetation file for Land Use.', &
2442	& "No land point in the map for point", &
2443	& '(verify your land use file.)')
2444	veget_next(ib,:) = veget_last(ib,:)
2445	ENDIF
2446
2447	ENDDO
2448	ENDIF
2449
2450	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_readvegetmax after updating'
2451	!
2452	frac_nobio_next (:,:) = un
2453	!
2454	!MM
2455	! Work only for one nnobio !! (ie ice)
2456	DO inobio=1,nnobio
2457	DO jv=1,nvm
2458	DO ib = 1, nbpt
2459	frac_nobio_next(ib,inobio) = frac_nobio_next(ib,inobio) - veget_next(ib,jv)
2460	ENDDO
2461	ENDDO
2462	ENDDO
2463
2464	DO ib = 1, nbpt
2465	sum_veg = SUM(veget_next(ib,:))
2466	sum_nobio = SUM(frac_nobio_next(ib,:))
2467	IF (sum_nobio < 0.) THEN
2468	frac_nobio_next(ib,:) = zero
2469	veget_next(ib,1) = veget_next(ib,1) + sum_nobio
2470	sum_veg = SUM(veget_next(ib,:))
2471	ENDIF
2472	sumf = sum_veg + sum_nobio
2473	IF (sumf > min_sechiba) THEN
2474	veget_next(ib,:) = veget_next(ib,:) / sumf
2475	frac_nobio_next(ib,:) = frac_nobio_next(ib,:) / sumf
2476	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2477	err=norm-un
2478	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ",ib, &
2479	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un, sumf",err,sumf
2480	IF (abs(err) > -EPSILON(un)) THEN
2481	IF ( SUM(frac_nobio_next(ib,:)) > min_sechiba ) THEN
2482	frac_nobio_next(ib,1) = frac_nobio_next(ib,1) - err
2483	ELSE
2484	veget_next(ib,1) = veget_next(ib,1) - err
2485	ENDIF
2486	norm=SUM(veget_next(ib,:))+SUM(frac_nobio_next(ib,:))
2487	err=norm-un
2488	IF (printlev_loc >=5) WRITE(numout,*) " slowproc_readvegetmax: ib ", ib, &
2489	" SUM(veget_next(ib,:)+frac_nobio_next(ib,:))-un",err
2490	IF (abs(err) > EPSILON(un)) THEN
2491	WRITE(numout,*) ' slowproc_readvegetmax _______'
2492	WRITE(numout,*) "update : Problem with point ",ib,",(",lalo(ib,1),",",lalo(ib,2),")"
2493	WRITE(numout,*) " err(sum-1.) = ",abs(err)
2494	CALL ipslerr_p (2,'slowproc_readvegetmax', &
2495	& 'Problem with sum vegetation + sum fracnobio for Land Use.', &
2496	& "sum not equal to 1.", &
2497	& '(verify your land use file.)')
2498	aveget(ib) = -0.6
2499	ENDIF
2500	ENDIF
2501	ELSE
2502	! sumf < min_sechiba
2503	WRITE(numout,*) ' slowproc_readvegetmax _______'
2504	WRITE(numout,*)" No vegetation nor frac_nobio for point ", ib, ",(", lalo(ib,1), ",", &
2505	lalo(ib,2),")"
2506	WRITE(numout,*)" Replaced by bare_soil !! "
2507	veget_next(ib,1) = un
2508	veget_next(ib,2:nvm) = zero
2509	frac_nobio_next(ib,:) = zero
2510	!!!$ CALL ipslerr_p (3,'slowproc_readvegetmax', &
2511	!!!$ & 'Problem with vegetation file for Land Use.', &
2512	!!!$ & "No vegetation nor frac_nobio for point ", &
2513	!!!$ & '(verify your land use file.)')
2514	ENDIF
2515	ENDDO
2516
2517	!! Set to zero fractions of frac_nobio and veget_max smaller than min_vegfrac
2518	!! Normalize to have the sum equal 1.
2519	CALL slowproc_veget_max_limit(nbpt, frac_nobio_next, veget_next)
2520
2521	! Write diagnostics
2522	CALL xios_orchidee_send_field("aveget",aveget)
2523	CALL xios_orchidee_send_field("interp_diag_aveget",aveget)
2524	CALL xios_orchidee_send_field("interp_diag_vegetrefrac",vegetrefrac)
2525
2526	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_readvegetmax ended'
2527
2528	END SUBROUTINE slowproc_readvegetmax
2529
2530
2531	!! ================================================================================================================================
2532	!! SUBROUTINE : slowproc_nearest
2533	!!
2534	!>\BRIEF looks for nearest grid point on the fine map
2535	!!
2536	!! DESCRIPTION : (definitions, functional, design, flags):
2537	!!
2538	!! RECENT CHANGE(S): None
2539	!!
2540	!! MAIN OUTPUT VARIABLE(S): ::inear
2541	!!
2542	!! REFERENCE(S) : None
2543	!!
2544	!! FLOWCHART : None
2545	!! \n
2546	!_ ================================================================================================================================
2547
2548	SUBROUTINE slowproc_nearest(iml, lon5, lat5, lonmod, latmod, inear)
2549
2550	!! INTERFACE DESCRIPTION
2551
2552	!! 0.1 input variables
2553
2554	INTEGER(i_std), INTENT(in) :: iml !! size of the vector
2555	REAL(r_std), DIMENSION(iml), INTENT(in) :: lon5, lat5 !! longitude and latitude vector, for the 5km vegmap
2556	REAL(r_std), INTENT(in) :: lonmod, latmod !! longitude and latitude modelled
2557
2558	!! 0.2 output variables
2559
2560	INTEGER(i_std), INTENT(out) :: inear !! location of the grid point from the 5km vegmap grid
2561	!! closest from the modelled grid point
2562
2563	!! 0.4 Local variables
2564
2565	REAL(r_std) :: pa, p
2566	REAL(r_std) :: coscolat, sincolat
2567	REAL(r_std) :: cospa, sinpa
2568	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: cosang
2569	INTEGER(i_std) :: i
2570	INTEGER(i_std), DIMENSION(1) :: ineartab
2571	INTEGER :: ALLOC_ERR
2572
2573	!_ ================================================================================================================================
2574
2575	ALLOCATE(cosang(iml), STAT=ALLOC_ERR)
2576	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_nearest','Error in allocation for cosang','','')
2577
2578	pa = pi/2.0 - latmod*pi/180.0 ! dist. between north pole and the point a
2579	!! COLATITUDE, in radian
2580	cospa = COS(pa)
2581	sinpa = SIN(pa)
2582
2583	DO i = 1, iml
2584
2585	sincolat = SIN( pi/2.0 - lat5(i)*pi/180.0 ) !! sinus of the colatitude
2586	coscolat = COS( pi/2.0 - lat5(i)*pi/180.0 ) !! cosinus of the colatitude
2587
2588	p = (lonmod-lon5(i))*pi/180.0 !! angle between a & b (between their meridian)in radians
2589
2590	!! dist(i) = ACOS( cospacoscolat + sinpasincolat*COS(p))
2591	cosang(i) = cospacoscolat + sinpasincolat*COS(p) !! TL : cosang is maximum when angle is at minimal value
2592	!! orthodromic distance between 2 points : cosang = cosinus (arc(AB)/R), with
2593	!R = Earth radius, then max(cosang) = max(cos(arc(AB)/R)), reached when arc(AB)/R is minimal, when
2594	! arc(AB) is minimal, thus when point B (corresponding grid point from LAI MAP) is the nearest from
2595	! modelled A point
2596	ENDDO
2597
2598	ineartab = MAXLOC( cosang(:) )
2599	inear = ineartab(1)
2600
2601	DEALLOCATE(cosang)
2602	END SUBROUTINE slowproc_nearest
2603
2604	!! ================================================================================================================================
2605	!! SUBROUTINE : slowproc_soilt
2606	!!
2607	!>\BRIEF Interpolate the Zobler or Reynolds/USDA soil type map
2608	!!
2609	!! DESCRIPTION : (definitions, functional, design, flags):
2610	!!
2611	!! RECENT CHANGE(S): Nov 2014, ADucharne
2612	!!
2613	!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction
2614	!!
2615	!! REFERENCE(S) : Reynold, Jackson, and Rawls (2000). Estimating soil water-holding capacities
2616	!! by linking the Food and Agriculture Organization soil map of the world with global pedon
2617	!! databases and continuous pedotransfer functions, WRR, 36, 3653-3662
2618	!!
2619	!! FLOWCHART : None
2620	!! \n
2621	!_ ================================================================================================================================
2622	SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction)
2623
2624	USE interpweight
2625
2626	IMPLICIT NONE
2627	!
2628	!
2629	! This subroutine should read the Zobler/Reynolds map and interpolate to the model grid.
2630	! The method is to get fraction of the three/12 main soiltypes for each grid box.
2631	! For the Zobler case, also called FAO in the code, the soil fraction are going to be put
2632	! into the array soiltype in the following order : coarse, medium and fine.
2633	! For the Reynolds/USDA case, the soiltype array follows the order defined in constantes_soil_var.f90
2634	!
2635	!
2636	!! 0.1 INPUT
2637	!
2638	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
2639	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
2640	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
2641	!! (1=North and then clockwise)
2642	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
2643	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
2644	!
2645	! 0.2 OUTPUT
2646	!
2647	REAL(r_std), INTENT(out) :: soilclass(nbpt, nscm) !! Soil type map to be created from the Zobler map
2648	!! or a map defining the 12 USDA classes (e.g. Reynolds)
2649	!! Holds the area of each texture class in the ORCHIDEE grid cells
2650	!! Final unit = fraction of ORCHIDEE grid-cell (unitless)
2651	REAL(r_std), INTENT(out) :: clayfraction(nbpt) !! The fraction of clay as used by STOMATE
2652	REAL(r_std), INTENT(out) :: sandfraction(nbpt) !! The fraction of sand (for SP-MIP)
2653	REAL(r_std), INTENT(out) :: siltfraction(nbpt) !! The fraction of silt (for SP-MIP)
2654	!
2655	!
2656	! 0.3 LOCAL
2657	!
2658	CHARACTER(LEN=80) :: filename
2659	INTEGER(i_std) :: ib, ilf, nbexp, i
2660	INTEGER(i_std) :: fopt !! Nb of pts from the texture map within one ORCHIDEE grid-cell
2661	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: solt !! Texture the different points from the input texture map
2662	!! in one ORCHIDEE grid cell (unitless)
2663	!
2664	! Number of texture classes in Zobler
2665	!
2666	INTEGER(i_std), PARAMETER :: nzobler = 7 !! Nb of texture classes according in the Zobler map
2667	REAL(r_std),ALLOCATABLE :: textfrac_table(:,:) !! conversion table between the texture index
2668	!! and the granulometric composition
2669	!
2670	INTEGER :: ALLOC_ERR
2671	INTEGER :: ntextinfile !! number of soil textures in the in the file
2672	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: textrefrac !! text fractions re-dimensioned
2673	REAL(r_std), DIMENSION(nbpt) :: atext !! Availability of the texture interpolation
2674	REAL(r_std) :: vmin, vmax !! min/max values to use for the
2675
2676	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
2677	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in input file
2678	REAL(r_std), DIMENSION(:), ALLOCATABLE :: variabletypevals !! Values for all the types of the variable
2679	!! (variabletypevals(1) = -un, not used)
2680	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
2681	!! 'XYKindTime': Input values are kinds
2682	!! of something with a temporal
2683	!! evolution on the dx*dy matrix'
2684	LOGICAL :: nonegative !! whether negative values should be removed
2685	CHARACTER(LEN=50) :: maskingtype !! Type of masking
2686	!! 'nomask': no-mask is applied
2687	!! 'mbelow': take values below maskvals(1)
2688	!! 'mabove': take values above maskvals(1)
2689	!! 'msumrange': take values within 2 ranges;
2690	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
2691	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
2692	!! (normalized by maskvals(3))
2693	!! 'var': mask values are taken from a
2694	!! variable inside the file (>0)
2695	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
2696	!! `maskingtype')
2697	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
2698	INTEGER(i_std), DIMENSION(:), ALLOCATABLE :: vecpos
2699	REAL(r_std) :: sgn !! sum of fractions excluding glaciers and ocean
2700	!_ ================================================================================================================================
2701
2702	IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt'
2703	!
2704	! Needs to be a configurable variable
2705	!
2706	!
2707	!Config Key = SOILCLASS_FILE
2708	!Config Desc = Name of file from which soil types are read
2709	!Config Def = soils_param.nc
2710	!Config If = NOT(IMPOSE_VEG)
2711	!Config Help = The name of the file to be opened to read the soil types.
2712	!Config The data from this file is then interpolated to the grid of
2713	!Config of the model. The aim is to get fractions for sand loam and
2714	!Config clay in each grid box. This information is used for soil hydrology
2715	!Config and respiration.
2716	!Config Units = [FILE]
2717	!
2718	! soils_param.nc file is 1deg soil texture file (Zobler)
2719	! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution)
2720
2721	filename = 'soils_param.nc'
2722	CALL getin_p('SOILCLASS_FILE',filename)
2723
2724	variablename = 'soiltext'
2725
2726	!! Variables for interpweight
2727	! Type of calculation of cell fractions
2728	fractype = 'default'
2729
2730	IF (xios_interpolation) THEN
2731	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " &
2732	// TRIM(filename) // " for variable " // TRIM(variablename)
2733
2734	SELECT CASE(soil_classif)
2735
2736	CASE('none')
2737	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
2738	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
2739	DO ib=1, nbpt
2740	soilclass(ib,:) = soilclass_default_fao
2741	clayfraction(ib) = clayfraction_default
2742	ENDDO
2743
2744
2745	CASE('zobler')
2746
2747	!
2748	soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
2749	!
2750	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
2751	!
2752	ALLOCATE(textrefrac(nbpt,nzobler))
2753	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
2754	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
2755	CALL get_soilcorr_zobler (nzobler, textfrac_table)
2756
2757	CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
2758	CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
2759	CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
2760	CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
2761	CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
2762	CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
2763	CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
2764
2765
2766
2767	CALL get_soilcorr_zobler (nzobler, textfrac_table)
2768	!
2769	!
2770	DO ib =1, nbpt
2771	soilclass(ib,1)=textrefrac(ib,1)
2772	soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
2773	soilclass(ib,3)=textrefrac(ib,5)
2774
2775	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
2776	! over the zobler pixels composing the ORCHIDEE grid-cell
2777	clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + &
2778	textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + &
2779	textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7)
2780
2781	sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + &
2782	textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + &
2783	textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7)
2784
2785	siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + &
2786	textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + &
2787	textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7)
2788
2789	sgn=SUM(soilclass(ib,1:3))
2790
2791	IF (sgn < min_sechiba) THEN
2792	soilclass(ib,:) = soilclass_default(:)
2793	clayfraction(ib) = clayfraction_default
2794	sandfraction(ib) = sandfraction_default
2795	siltfraction(ib) = siltfraction_default
2796	atext(ib)=0.
2797	ELSE
2798	atext(ib)=sgn
2799	clayfraction(ib) = clayfraction(ib) / sgn
2800	sandfraction(ib) = sandfraction(ib) / sgn
2801	siltfraction(ib) = siltfraction(ib) / sgn
2802	soilclass(ib,1:3) = soilclass(ib,1:3) / sgn
2803	ENDIF
2804
2805	ENDDO
2806
2807
2808
2809	CASE('usda')
2810
2811	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda'
2812
2813	soilclass_default=soilclass_default_usda
2814	!
2815	WRITE(numout,*) "Using a soilclass map with usda classification"
2816	!
2817	ALLOCATE(textrefrac(nbpt,nscm))
2818	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
2819	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
2820
2821	CALL get_soilcorr_usda (nscm, textfrac_table)
2822
2823	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
2824
2825	CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
2826	CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
2827	CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
2828	CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
2829	CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
2830	CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
2831	CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
2832	CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8))
2833	CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9))
2834	CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10))
2835	CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11))
2836	CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12))
2837
2838
2839
2840	CALL get_soilcorr_usda (nscm, textfrac_table)
2841	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
2842
2843	DO ib =1, nbpt
2844	clayfraction(ib) = 0.0
2845	DO ilf = 1,nscm
2846	soilclass(ib,ilf)=textrefrac(ib,ilf)
2847	clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf)
2848	sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf)
2849	siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf)
2850	ENDDO
2851
2852
2853	sgn=SUM(soilclass(ib,:))
2854
2855	IF (sgn < min_sechiba) THEN
2856	soilclass(ib,:) = soilclass_default(:)
2857	clayfraction(ib) = clayfraction_default
2858	sandfraction(ib) = sandfraction_default
2859	siltfraction(ib) = siltfraction_default
2860	atext(ib)=0
2861	ELSE
2862	soilclass(ib,:) = soilclass(ib,:) / sgn
2863	clayfraction(ib) = clayfraction(ib) / sgn
2864	sandfraction(ib) = sandfraction(ib) / sgn
2865	siltfraction(ib) = siltfraction(ib) / sgn
2866	atext(ib)=sgn
2867	ENDIF
2868	ENDDO
2869
2870	CASE DEFAULT
2871	WRITE(numout,*) 'slowproc_soilt:'
2872	WRITE(numout,*) ' A non supported soil type classification has been chosen'
2873	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
2874	END SELECT
2875
2876
2877
2878	ELSE ! xios_interpolation
2879	! Read and interpolate using stardard method with IOIPSL and aggregate
2880
2881	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
2882	// TRIM(filename) // " for variable " // TRIM(variablename)
2883
2884
2885	! Name of the longitude and latitude in the input file
2886	lonname = 'nav_lon'
2887	latname = 'nav_lat'
2888
2889	IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " &
2890	// TRIM(filename) // " for variable " // TRIM(variablename)
2891
2892	IF ( TRIM(soil_classif) /= 'none' ) THEN
2893
2894	! Define a variable for the number of soil textures in the input file
2895	SELECTCASE(soil_classif)
2896	CASE('zobler')
2897	ntextinfile=nzobler
2898	CASE('usda')
2899	ntextinfile=nscm
2900	CASE DEFAULT
2901	WRITE(numout,*) 'slowproc_soilt:'
2902	WRITE(numout,*) ' A non supported soil type classification has been chosen'
2903	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
2904	ENDSELECT
2905
2906	ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR)
2907	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',&
2908	'','')
2909
2910	! Assigning values to vmin, vmax
2911	vmin = un
2912	vmax = ntextinfile*un
2913
2914	ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR)
2915	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','')
2916	variabletypevals = -un
2917
2918	!! Variables for interpweight
2919	! Should negative values be set to zero from input file?
2920	nonegative = .FALSE.
2921	! Type of mask to apply to the input data (see header for more details)
2922	maskingtype = 'mabove'
2923	! Values to use for the masking
2924	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
2925	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
2926	namemaskvar = ''
2927
2928	CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours, &
2929	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
2930	maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext)
2931
2932	ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR)
2933	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','')
2934	ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR)
2935	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','')
2936
2937	IF (printlev_loc >= 5) THEN
2938	WRITE(numout,*)' slowproc_soilt after interpweight_2D'
2939	WRITE(numout,*)' slowproc_soilt before starting loop nbpt:', nbpt
2940	WRITE(numout,*)" slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..."
2941	END IF
2942	ELSE
2943	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt using default values all points are propertly ' // &
2944	'interpolated atext = 1. everywhere!'
2945	atext = 1.
2946	END IF
2947
2948	nbexp = 0
2949	SELECTCASE(soil_classif)
2950	CASE('none')
2951	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
2952	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
2953	DO ib=1, nbpt
2954	soilclass(ib,:) = soilclass_default_fao
2955	clayfraction(ib) = clayfraction_default
2956	sandfraction(ib) = sandfraction_default
2957	siltfraction(ib) = siltfraction_default
2958	ENDDO
2959	CASE('zobler')
2960	!
2961	soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
2962	!
2963	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
2964	!
2965	ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
2966	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
2967	CALL get_soilcorr_zobler (nzobler, textfrac_table)
2968	!
2969	!
2970	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_soilt after getting table of textures'
2971	DO ib =1, nbpt
2972	soilclass(ib,:) = zero
2973	clayfraction(ib) = zero
2974	sandfraction(ib) = zero
2975	siltfraction(ib) = zero
2976	!
2977	! vecpos: List of positions where textures were not zero
2978	! vecpos(1): number of not null textures found
2979	vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq')
2980	fopt = vecpos(1)
2981
2982	IF ( fopt .EQ. 0 ) THEN
2983	! No points were found for current grid box, use default values
2984	nbexp = nbexp + 1
2985	soilclass(ib,:) = soilclass_default(:)
2986	clayfraction(ib) = clayfraction_default
2987	sandfraction(ib) = sandfraction_default
2988	siltfraction(ib) = siltfraction_default
2989	ELSE
2990	IF (fopt == nzobler) THEN
2991	! All textures are not zero
2992	solt=(/(i,i=1,nzobler)/)
2993	ELSE
2994	DO ilf = 1,fopt
2995	solt(ilf) = vecpos(ilf+1)
2996	END DO
2997	END IF
2998	!
2999	! Compute the fraction of each textural class
3000	!
3001	sgn = 0.
3002	DO ilf = 1,fopt
3003	!
3004	! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE
3005	! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine)
3006	! For type 6 = glacier, default values are set and it is also taken into account during the normalization
3007	! of the fractions (done in interpweight_2D)
3008	! Note that type 0 corresponds to ocean but it is already removed using the mask above.
3009	!
3010	IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. &
3011	(solt(ilf) .NE. 6) ) THEN
3012	SELECT CASE(solt(ilf))
3013	CASE(1)
3014	soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf))
3015	CASE(2)
3016	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3017	CASE(3)
3018	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3019	CASE(4)
3020	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3021	CASE(5)
3022	soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf))
3023	CASE(7)
3024	soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
3025	CASE DEFAULT
3026	WRITE(numout,*) 'We should not be here, an impossible case appeared'
3027	CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','')
3028	END SELECT
3029	! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
3030	! over the zobler pixels composing the ORCHIDEE grid-cell
3031	clayfraction(ib) = clayfraction(ib) + &
3032	& textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf))
3033	sandfraction(ib) = sandfraction(ib) + &
3034	& textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf))
3035	siltfraction(ib) = siltfraction(ib) + &
3036	& textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf))
3037	! Sum the fractions which are not glaciers nor ocean
3038	sgn = sgn + textrefrac(ib,solt(ilf))
3039	ELSE
3040	IF (solt(ilf) .GT. nzobler) THEN
3041	WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
3042	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','')
3043	ENDIF
3044	END IF
3045	ENDDO
3046
3047	IF ( sgn .LT. min_sechiba) THEN
3048	! Set default values if grid cells were only covered by glaciers or ocean
3049	! or if now information on the source grid was found.
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	! Normalize using the fraction of surface not including glaciers and ocean
3057	soilclass(ib,:) = soilclass(ib,:)/sgn
3058	clayfraction(ib) = clayfraction(ib)/sgn
3059	sandfraction(ib) = sandfraction(ib)/sgn
3060	siltfraction(ib) = siltfraction(ib)/sgn
3061	ENDIF
3062	ENDIF
3063	ENDDO
3064
3065	! The "USDA" case reads a map of the 12 USDA texture classes,
3066	! such as to assign the corresponding soil properties
3067	CASE("usda")
3068	IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification"
3069
3070	soilclass_default=soilclass_default_usda
3071
3072	ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
3073	IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
3074
3075	CALL get_soilcorr_usda (nscm, textfrac_table)
3076
3077	IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
3078	!
3079	DO ib =1, nbpt
3080	! GO through the point we have found
3081	!
3082	!
3083	! Provide which textures were found
3084	! vecpos: List of positions where textures were not zero
3085	! vecpos(1): number of not null textures found
3086	vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq')
3087	fopt = vecpos(1)
3088
3089	!
3090	! Check that we found some points
3091	!
3092	soilclass(ib,:) = 0.0
3093	clayfraction(ib) = 0.0
3094
3095	IF ( fopt .EQ. 0) THEN
3096	! No points were found for current grid box, use default values
3097	IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib
3098	nbexp = nbexp + 1
3099	soilclass(ib,:) = soilclass_default
3100	clayfraction(ib) = clayfraction_default
3101	sandfraction(ib) = sandfraction_default
3102	siltfraction(ib) = siltfraction_default
3103	ELSE
3104	IF (fopt == nscm) THEN
3105	! All textures are not zero
3106	solt(:) = (/(i,i=1,nscm)/)
3107	ELSE
3108	DO ilf = 1,fopt
3109	solt(ilf) = vecpos(ilf+1)
3110	END DO
3111	END IF
3112
3113	!
3114	!
3115	! Compute the fraction of each textural class
3116	!
3117	!
3118	DO ilf = 1,fopt
3119	IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
3120	soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf))
3121	clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * &
3122	textrefrac(ib,solt(ilf))
3123	sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * &
3124	textrefrac(ib,solt(ilf))
3125	siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * &
3126	textrefrac(ib,solt(ilf))
3127	ELSE
3128	IF (solt(ilf) .GT. nscm) THEN
3129	WRITE(,) 'The file contains a soil color class which is incompatible with this program'
3130	CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','')
3131	ENDIF
3132	ENDIF
3133	!
3134	ENDDO
3135
3136	! Set default values if the surface in source file is too small
3137	IF ( atext(ib) .LT. min_sechiba) THEN
3138	nbexp = nbexp + 1
3139	soilclass(ib,:) = soilclass_default(:)
3140	clayfraction(ib) = clayfraction_default
3141	sandfraction(ib) = sandfraction_default
3142	siltfraction(ib) = siltfraction_default
3143	ENDIF
3144	ENDIF
3145
3146	ENDDO
3147
3148	IF (printlev_loc>=4) WRITE (numout,*) ' slowproc_soilt: End case usda'
3149
3150	CASE DEFAULT
3151	WRITE(numout,*) 'slowproc_soilt _______'
3152	WRITE(numout,*) ' A non supported soil type classification has been chosen'
3153	CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
3154	ENDSELECT
3155	IF (printlev_loc >= 5 ) WRITE(numout,*)' slowproc_soilt end of type classification'
3156
3157	IF ( nbexp .GT. 0 ) THEN
3158	WRITE(numout,*) 'slowproc_soilt:'
3159	WRITE(numout,*) ' The interpolation of the bare soil albedo had ', nbexp
3160	WRITE(numout,*) ' points without data. This are either coastal points or ice covered land.'
3161	WRITE(numout,*) ' The problem was solved by using the default soil types.'
3162	ENDIF
3163
3164	IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals)
3165	IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac)
3166	IF (ALLOCATED(solt)) DEALLOCATE (solt)
3167	IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table)
3168
3169	ENDIF ! xios_interpolation
3170
3171	! Write diagnostics
3172	CALL xios_orchidee_send_field("atext",atext)
3173
3174	CALL xios_orchidee_send_field("interp_diag_atext",atext)
3175	CALL xios_orchidee_send_field("interp_diag_soilclass",soilclass)
3176	CALL xios_orchidee_send_field("interp_diag_clayfraction",clayfraction)
3177
3178	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_soilt ended'
3179
3180	END SUBROUTINE slowproc_soilt
3181
3182	!! ================================================================================================================================
3183	!! SUBROUTINE : slowproc_slope
3184	!!
3185	!>\BRIEF Calculate mean slope coef in each model grid box from the slope map
3186	!!
3187	!! DESCRIPTION : (definitions, functional, design, flags):
3188	!!
3189	!! RECENT CHANGE(S): None
3190	!!
3191	!! MAIN OUTPUT VARIABLE(S): ::reinf_slope
3192	!!
3193	!! REFERENCE(S) : None
3194	!!
3195	!! FLOWCHART : None
3196	!! \n
3197	!_ ================================================================================================================================
3198
3199	SUBROUTINE slowproc_slope(nbpt, lalo, neighbours, resolution, contfrac, reinf_slope)
3200
3201	USE interpweight
3202
3203	IMPLICIT NONE
3204
3205	!
3206	!
3207	!
3208	! 0.1 INPUT
3209	!
3210	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3211	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
3212	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point
3213	! (1=North and then clockwise)
3214	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
3215	REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Fraction of continent in the grid
3216	!
3217	! 0.2 OUTPUT
3218	!
3219	REAL(r_std), INTENT(out) :: reinf_slope(nbpt) ! slope coef
3220	!
3221	! 0.3 LOCAL
3222	!
3223	!
3224	REAL(r_std) :: slope_noreinf ! Slope above which runoff is maximum
3225	CHARACTER(LEN=80) :: filename
3226	REAL(r_std) :: vmin, vmax !! min/max values to use for the
3227	!! renormalization
3228	REAL(r_std), DIMENSION(nbpt) :: aslope !! slope availability
3229
3230	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3231	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
3232	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3233	!! 'XYKindTime': Input values are kinds
3234	!! of something with a temporal
3235	!! evolution on the dx*dy matrix'
3236	LOGICAL :: nonegative !! whether negative values should be removed
3237	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3238	!! 'nomask': no-mask is applied
3239	!! 'mbelow': take values below maskvals(1)
3240	!! 'mabove': take values above maskvals(1)
3241	!! 'msumrange': take values within 2 ranges;
3242	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3243	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3244	!! (normalized by maskvals(3))
3245	!! 'var': mask values are taken from a
3246	!! variable inside the file (>0)
3247	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3248	!! `maskingtype')
3249	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3250
3251	!_ ================================================================================================================================
3252
3253	!
3254	!Config Key = SLOPE_NOREINF
3255	!Config Desc = See slope_noreinf above
3256	!Config If =
3257	!Config Def = 0.5
3258	!Config Help = The slope above which there is no reinfiltration
3259	!Config Units = [-]
3260	!
3261	slope_noreinf = 0.5
3262	!
3263	CALL getin_p('SLOPE_NOREINF',slope_noreinf)
3264	!
3265	!Config Key = TOPOGRAPHY_FILE
3266	!Config Desc = Name of file from which the topography map is to be read
3267	!Config If =
3268	!Config Def = cartepente2d_15min.nc
3269	!Config Help = The name of the file to be opened to read the orography
3270	!Config map is to be given here. Usualy SECHIBA runs with a 2'
3271	!Config map which is derived from the NGDC one.
3272	!Config Units = [FILE]
3273	!
3274	filename = 'cartepente2d_15min.nc'
3275	CALL getin_p('TOPOGRAPHY_FILE',filename)
3276
3277	IF (xios_interpolation) THEN
3278
3279	CALL xios_orchidee_recv_field('reinf_slope_interp',reinf_slope)
3280	CALL xios_orchidee_recv_field('frac_slope_interp',aslope)
3281
3282
3283	ELSE
3284
3285	variablename = 'pente'
3286	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_slope: Read and interpolate " &
3287	// TRIM(filename) // " for variable " // TRIM(variablename)
3288
3289	! For this case there are not types/categories. We have 'only' a continuos field
3290	! Assigning values to vmin, vmax
3291	vmin = 0.
3292	vmax = 9999.
3293
3294	!! Variables for interpweight
3295	! Type of calculation of cell fractions
3296	fractype = 'slopecalc'
3297	! Name of the longitude and latitude in the input file
3298	lonname = 'longitude'
3299	latname = 'latitude'
3300	! Should negative values be set to zero from input file?
3301	nonegative = .FALSE.
3302	! Type of mask to apply to the input data (see header for more details)
3303	maskingtype = 'mabove'
3304	! Values to use for the masking
3305	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3306	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
3307	namemaskvar = ''
3308
3309	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3310	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3311	maskvals, namemaskvar, -1, fractype, slope_default, slope_noreinf, &
3312	reinf_slope, aslope)
3313	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_slope after interpweight_2Dcont'
3314
3315	ENDIF
3316
3317	! Write diagnostics
3318	CALL xios_orchidee_send_field("aslope",aslope)
3319	CALL xios_orchidee_send_field("interp_diag_aslope",aslope)
3320
3321	CALL xios_orchidee_send_field("interp_diag_reinf_slope",reinf_slope)
3322
3323	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_slope ended'
3324
3325	END SUBROUTINE slowproc_slope
3326
3327
3328	!! ================================================================================================================================
3329	!! SUBROUTINE : slowproc_woodharvest
3330	!!
3331	!>\BRIEF
3332	!!
3333	!! DESCRIPTION :
3334	!!
3335	!! RECENT CHANGE(S): None
3336	!!
3337	!! MAIN OUTPUT VARIABLE(S): ::
3338	!!
3339	!! REFERENCE(S) : None
3340	!!
3341	!! FLOWCHART : None
3342	!! \n
3343	!_ ================================================================================================================================
3344
3345	SUBROUTINE slowproc_woodharvest(nbpt, lalo, neighbours, resolution, contfrac, woodharvest)
3346
3347	USE interpweight
3348
3349	IMPLICIT NONE
3350
3351	!
3352	!
3353	!
3354	! 0.1 INPUT
3355	!
3356	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data needs to be interpolated
3357	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
3358	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
3359	!! (1=North and then clockwise)
3360	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
3361	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of continent in the grid
3362	!
3363	! 0.2 OUTPUT
3364	!
3365	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: woodharvest !! Wood harvest
3366	!
3367	! 0.3 LOCAL
3368	!
3369	CHARACTER(LEN=80) :: filename
3370	REAL(r_std) :: vmin, vmax
3371	REAL(r_std), DIMENSION(nbpt) :: aoutvar !! availability of input data to
3372	!! interpolate output variable
3373	!! (on the nbpt space)
3374	CHARACTER(LEN=80) :: variablename !! Variable to interpolate
3375	CHARACTER(LEN=80) :: lonname, latname !! lon, lat name in the input file
3376	CHARACTER(LEN=50) :: fractype !! method of calculation of fraction
3377	!! 'XYKindTime': Input values are kinds
3378	!! of something with a temporal
3379	!! evolution on the dx*dy matrix'
3380	LOGICAL :: nonegative !! whether negative values should be removed
3381	CHARACTER(LEN=50) :: maskingtype !! Type of masking
3382	!! 'nomask': no-mask is applied
3383	!! 'mbelow': take values below maskvals(1)
3384	!! 'mabove': take values above maskvals(1)
3385	!! 'msumrange': take values within 2 ranges;
3386	!! maskvals(2) <= SUM(vals(k)) <= maskvals(1)
3387	!! maskvals(1) < SUM(vals(k)) <= maskvals(3)
3388	!! (normalized by maskvals(3))
3389	!! 'var': mask values are taken from a
3390	!! variable inside the file (>0)
3391	REAL(r_std), DIMENSION(3) :: maskvals !! values to use to mask (according to
3392	!! `maskingtype')
3393	CHARACTER(LEN=250) :: namemaskvar !! name of the variable to use to mask
3394	REAL(r_std), DIMENSION(1) :: variabletypevals !!
3395	! REAL(r_std), DIMENSION(nbp_mpi) :: woodharvest_mpi !! Wood harvest where all thredds OMP are gatherd
3396	!_ ================================================================================================================================
3397
3398
3399	!Config Key = WOODHARVEST_FILE
3400	!Config Desc = Name of file from which the wood harvest will be read
3401	!Config If = DO_WOOD_HARVEST
3402	!Config Def = woodharvest.nc
3403	!Config Help =
3404	!Config Units = [FILE]
3405	filename = 'woodharvest.nc'
3406	CALL getin_p('WOODHARVEST_FILE',filename)
3407	variablename = 'woodharvest'
3408
3409
3410	IF (xios_interpolation) THEN
3411	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readwoodharvest: Use XIOS to read and interpolate " &
3412	// TRIM(filename) // " for variable " // TRIM(variablename)
3413
3414	CALL xios_orchidee_recv_field('woodharvest_interp',woodharvest)
3415
3416	aoutvar = 1.0
3417	ELSE
3418
3419	IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_readwoodharvest: Read and interpolate " &
3420	// TRIM(filename) // " for variable " // TRIM(variablename)
3421
3422	! For this case there are not types/categories. We have 'only' a continuos field
3423	! Assigning values to vmin, vmax
3424	vmin = 0.
3425	vmax = 9999.
3426
3427	!! Variables for interpweight
3428	! Type of calculation of cell fractions
3429	fractype = 'default'
3430	! Name of the longitude and latitude in the input file
3431	lonname = 'longitude'
3432	latname = 'latitude'
3433	! Should negative values be set to zero from input file?
3434	nonegative = .TRUE.
3435	! Type of mask to apply to the input data (see header for more details)
3436	maskingtype = 'nomask'
3437	! Values to use for the masking
3438	maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
3439	! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') (here not used)
3440	namemaskvar = ''
3441
3442	variabletypevals=-un
3443	CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours, &
3444	contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype, &
3445	maskvals, namemaskvar, -1, fractype, 0., 0., woodharvest, aoutvar)
3446	IF (printlev_loc >= 5) WRITE(numout,*)' slowproc_wodharvest after interpweight_2Dcont'
3447
3448	IF (printlev_loc >= 3) WRITE(numout,*) ' slowproc_woodharvest ended'
3449	END IF
3450	END SUBROUTINE slowproc_woodharvest
3451
3452
3453	!! ================================================================================================================================
3454	!! SUBROUTINE : get_soilcorr_zobler
3455	!!
3456	!>\BRIEF The "get_soilcorr" routine defines the table of correspondence
3457	!! between the Zobler types and the three texture types known by SECHIBA and STOMATE :
3458	!! silt, sand and clay.
3459	!!
3460	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
3461	!! The data from this file is then interpolated to the grid of the model. \n
3462	!! The aim is to get fractions for sand loam and clay in each grid box.\n
3463	!! This information is used for soil hydrology and respiration.
3464	!!
3465	!!
3466	!! RECENT CHANGE(S): None
3467	!!
3468	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
3469	!!
3470	!! REFERENCE(S) :
3471	!! - Zobler L., 1986, A World Soil File for global climate modelling. NASA Technical memorandum 87802. NASA
3472	!! Goddard Institute for Space Studies, New York, U.S.A.
3473	!!
3474	!! FLOWCHART : None
3475	!! \n
3476	!_ ================================================================================================================================
3477
3478	SUBROUTINE get_soilcorr_zobler (nzobler,textfrac_table)
3479
3480	IMPLICIT NONE
3481
3482	!! 0. Variables and parameters declaration
3483
3484	INTEGER(i_std),PARAMETER :: nbtypes_zobler = 7 !! Number of Zobler types (unitless)
3485
3486	!! 0.1 Input variables
3487
3488	INTEGER(i_std),INTENT(in) :: nzobler !! Size of the array (unitless)
3489
3490	!! 0.2 Output variables
3491
3492	REAL(r_std),DIMENSION(nzobler,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
3493	!! and granulometric composition (0-1, unitless)
3494
3495	!! 0.4 Local variables
3496
3497	INTEGER(i_std) :: ib !! Indice (unitless)
3498
3499	!_ ================================================================================================================================
3500
3501	!-
3502	! 0. Check consistency
3503	!-
3504	IF (nzobler /= nbtypes_zobler) THEN
3505	CALL ipslerr_p(3,'get_soilcorr', 'nzobler /= nbtypes_zobler',&
3506	& 'We do not have the correct number of classes', &
3507	& ' in the code for the file.') ! Fatal error
3508	ENDIF
3509
3510	!-
3511	! 1. Textural fraction for : silt sand clay
3512	!-
3513	textfrac_table(1,:) = (/ 0.12, 0.82, 0.06 /)
3514	textfrac_table(2,:) = (/ 0.32, 0.58, 0.10 /)
3515	textfrac_table(3,:) = (/ 0.39, 0.43, 0.18 /)
3516	textfrac_table(4,:) = (/ 0.15, 0.58, 0.27 /)
3517	textfrac_table(5,:) = (/ 0.34, 0.32, 0.34 /)
3518	textfrac_table(6,:) = (/ 0.00, 1.00, 0.00 /)
3519	textfrac_table(7,:) = (/ 0.39, 0.43, 0.18 /)
3520
3521
3522	!-
3523	! 2. Check the mapping for the Zobler types which are going into the ORCHIDEE textures classes
3524	!-
3525	DO ib=1,nzobler ! Loop over # classes soil
3526
3527	IF (ABS(SUM(textfrac_table(ib,:))-1.0) > EPSILON(1.0)) THEN ! The sum of the textural fractions should not exceed 1 !
3528	WRITE(numout,*) &
3529	& 'Error in the correspondence table', &
3530	& ' sum is not equal to 1 in', ib
3531	WRITE(numout,*) textfrac_table(ib,:)
3532	CALL ipslerr_p(3,'get_soilcorr', 'SUM(textfrac_table(ib,:)) /= 1.0',&
3533	& '', 'Error in the correspondence table') ! Fatal error
3534	ENDIF
3535
3536	ENDDO ! Loop over # classes soil
3537
3538
3539	END SUBROUTINE get_soilcorr_zobler
3540
3541	!! ================================================================================================================================
3542	!! SUBROUTINE : get_soilcorr_usda
3543	!!
3544	!>\BRIEF The "get_soilcorr_usda" routine defines the table of correspondence
3545	!! between the 12 USDA textural classes and their granulometric composition,
3546	!! as % of silt, sand and clay. This is used to further defien clayfraction.
3547	!!
3548	!! DESCRIPTION : get_soilcorr is needed if you use soils_param.nc .\n
3549	!! The data from this file is then interpolated to the grid of the model. \n
3550	!! The aim is to get fractions for sand loam and clay in each grid box.\n
3551	!! This information is used for soil hydrology and respiration.
3552	!! The default map in this case is derived from Reynolds et al 2000, \n
3553	!! at the 1/12deg resolution, with indices that are consistent with the \n
3554	!! textures tabulated below
3555	!!
3556	!! RECENT CHANGE(S): Created by A. Ducharne on July 02, 2014
3557	!!
3558	!! MAIN OUTPUT VARIABLE(S) : ::texfrac_table
3559	!!
3560	!! REFERENCE(S) :
3561	!!
3562	!! FLOWCHART : None
3563	!! \n
3564	!_ ================================================================================================================================
3565
3566	SUBROUTINE get_soilcorr_usda (nusda,textfrac_table)
3567
3568	IMPLICIT NONE
3569
3570	!! 0. Variables and parameters declaration
3571
3572	!! 0.1 Input variables
3573
3574	INTEGER(i_std),INTENT(in) :: nusda !! Size of the array (unitless)
3575
3576	!! 0.2 Output variables
3577
3578	REAL(r_std),DIMENSION(nusda,ntext),INTENT(out) :: textfrac_table !! Table of correspondence between soil texture class
3579	!! and granulometric composition (0-1, unitless)
3580
3581	!! 0.4 Local variables
3582
3583	INTEGER(i_std),PARAMETER :: nbtypes_usda = 12 !! Number of USDA texture classes (unitless)
3584	INTEGER(i_std) :: n !! Index (unitless)
3585
3586	!_ ================================================================================================================================
3587
3588	!-
3589	! 0. Check consistency
3590	!-
3591	IF (nusda /= nbtypes_usda) THEN
3592	CALL ipslerr_p(3,'get_soilcorr', 'nusda /= nbtypes_usda',&
3593	& 'We do not have the correct number of classes', &
3594	& ' in the code for the file.') ! Fatal error
3595	ENDIF
3596
3597	!! Parameters for soil type distribution :
3598	!! Sand, Loamy Sand, Sandy Loam, Silt Loam, Silt, Loam, Sandy Clay Loam, Silty Clay Loam, Clay Loam, Sandy Clay, Silty Clay, Clay
3599	! The order comes from constantes_soil.f90
3600	! The corresponding granulometric composition comes from Carsel & Parrish, 1988
3601
3602	!-
3603	! 1. Textural fractions for : sand, clay
3604	!-
3605	textfrac_table(1,2:3) = (/ 0.93, 0.03 /) ! Sand
3606	textfrac_table(2,2:3) = (/ 0.81, 0.06 /) ! Loamy Sand
3607	textfrac_table(3,2:3) = (/ 0.63, 0.11 /) ! Sandy Loam
3608	textfrac_table(4,2:3) = (/ 0.17, 0.19 /) ! Silt Loam
3609	textfrac_table(5,2:3) = (/ 0.06, 0.10 /) ! Silt
3610	textfrac_table(6,2:3) = (/ 0.40, 0.20 /) ! Loam
3611	textfrac_table(7,2:3) = (/ 0.54, 0.27 /) ! Sandy Clay Loam
3612	textfrac_table(8,2:3) = (/ 0.08, 0.33 /) ! Silty Clay Loam
3613	textfrac_table(9,2:3) = (/ 0.30, 0.33 /) ! Clay Loam
3614	textfrac_table(10,2:3) = (/ 0.48, 0.41 /) ! Sandy Clay
3615	textfrac_table(11,2:3) = (/ 0.06, 0.46 /) ! Silty Clay
3616	textfrac_table(12,2:3) = (/ 0.15, 0.55 /) ! Clay
3617
3618	! Fraction of silt
3619
3620	DO n=1,nusda
3621	textfrac_table(n,1) = 1. - textfrac_table(n,2) - textfrac_table(n,3)
3622	END DO
3623
3624	END SUBROUTINE get_soilcorr_usda
3625
3626	!! ================================================================================================================================
3627	!! FUNCTION : tempfunc
3628	!!
3629	!>\BRIEF ! This function interpolates value between ztempmin and ztempmax
3630	!! used for lai detection.
3631	!!
3632	!! DESCRIPTION : This subroutine calculates a scalar between 0 and 1 with the following equation :\n
3633	!! \latexonly
3634	!! \input{constantes_veg_tempfunc.tex}
3635	!! \endlatexonly
3636	!!
3637	!! RECENT CHANGE(S): None
3638	!!
3639	!! RETURN VALUE : tempfunc_result
3640	!!
3641	!! REFERENCE(S) : None
3642	!!
3643	!! FLOWCHART : None
3644	!! \n
3645	!_ ================================================================================================================================
3646
3647	FUNCTION tempfunc (temp_in) RESULT (tempfunc_result)
3648
3649
3650	!! 0. Variables and parameters declaration
3651
3652	REAL(r_std),PARAMETER :: ztempmin=273._r_std !! Temperature for laimin (K)
3653	REAL(r_std),PARAMETER :: ztempmax=293._r_std !! Temperature for laimax (K)
3654	REAL(r_std) :: zfacteur !! Interpolation factor (K^{-2})
3655
3656	!! 0.1 Input variables
3657
3658	REAL(r_std),INTENT(in) :: temp_in !! Temperature (K)
3659
3660	!! 0.2 Result
3661
3662	REAL(r_std) :: tempfunc_result !! (unitless)
3663
3664	!_ ================================================================================================================================
3665
3666	!! 1. Define a coefficient
3667	zfacteur = un/(ztempmax-ztempmin)**2
3668
3669	!! 2. Computes tempfunc
3670	IF (temp_in > ztempmax) THEN
3671	tempfunc_result = un
3672	ELSEIF (temp_in < ztempmin) THEN
3673	tempfunc_result = zero
3674	ELSE
3675	tempfunc_result = un-zfacteur(ztempmax-temp_in)*2
3676	ENDIF !(temp_in > ztempmax)
3677
3678
3679	END FUNCTION tempfunc
3680
3681
3682	!! ================================================================================================================================
3683	!! SUBROUTINE : slowproc_checkveget
3684	!!
3685	!>\BRIEF To verify the consistency of the various fractions defined within the grid box after having been
3686	!! been updated by STOMATE or the standard procedures.
3687	!!
3688	!! DESCRIPTION : (definitions, functional, design, flags):
3689	!!
3690	!! RECENT CHANGE(S): None
3691	!!
3692	!! MAIN OUTPUT VARIABLE(S): :: none
3693	!!
3694	!! REFERENCE(S) : None
3695	!!
3696	!! FLOWCHART : None
3697	!! \n
3698	!_ ================================================================================================================================
3699	!
3700	SUBROUTINE slowproc_checkveget(nbpt, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
3701
3702	! 0.1 INPUT
3703	!
3704	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
3705	REAL(r_std),DIMENSION (nbpt,nnobio), INTENT(in) :: frac_nobio ! Fraction of ice,lakes,cities, ... (unitless)
3706	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget_max ! Maximum fraction of vegetation type including none biological fraction (unitless)
3707	REAL(r_std),DIMENSION (nbpt,nvm), INTENT(in) :: veget ! Vegetation fractions
3708	REAL(r_std),DIMENSION (nbpt), INTENT(in) :: tot_bare_soil ! Total evaporating bare soil fraction within the mesh
3709	REAL(r_std),DIMENSION (nbpt,nstm), INTENT(in) :: soiltile ! Fraction of soil tiles in the gridbox (unitless)
3710
3711	! 0.3 LOCAL
3712	!
3713	INTEGER(i_std) :: ji, jn, jv
3714	REAL(r_std) :: epsilocal !! A very small value
3715	REAL(r_std) :: totfrac
3716	CHARACTER(len=80) :: str1, str2
3717
3718	!_ ================================================================================================================================
3719
3720	!
3721	! There is some margin added as the computing errors might bring us above EPSILON(un)
3722	!
3723	epsilocal = EPSILON(un)*1000.
3724
3725	!! 1.0 Verify that none of the fractions are smaller than min_vegfrac, without beeing zero.
3726	!!
3727	DO ji=1,nbpt
3728	DO jn=1,nnobio
3729	IF ( frac_nobio(ji,jn) > epsilocal .AND. frac_nobio(ji,jn) < min_vegfrac ) THEN
3730	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
3731	WRITE(str2,'("The small value obtained is ", E14.4)') frac_nobio(ji,jn)
3732	CALL ipslerr_p (3,'slowproc_checkveget', &
3733	"frac_nobio is larger than zero but smaller than min_vegfrac.", str1, str2)
3734	ENDIF
3735	ENDDO
3736	END DO
3737
3738	IF (.NOT. ok_dgvm) THEN
3739	DO ji=1,nbpt
3740	DO jv=1,nvm
3741	IF ( veget_max(ji,jv) > epsilocal .AND. veget_max(ji,jv) < min_vegfrac ) THEN
3742	WRITE(str1,'("Occurs on grid box", I8," and nobio type ",I3 )') ji, jn
3743	WRITE(str2,'("The small value obtained is ", E14.4)') veget_max(ji,jv)
3744	CALL ipslerr_p (3,'slowproc_checkveget', &
3745	"veget_max is larger than zero but smaller than min_vegfrac.", str1, str2)
3746	ENDIF
3747	ENDDO
3748	ENDDO
3749	END IF
3750
3751	!! 2.0 verify that with all the fractions we cover the entire grid box
3752	!!
3753	DO ji=1,nbpt
3754	totfrac = zero
3755	DO jn=1,nnobio
3756	totfrac = totfrac + frac_nobio(ji,jn)
3757	ENDDO
3758	DO jv=1,nvm
3759	totfrac = totfrac + veget_max(ji,jv)
3760	ENDDO
3761	IF ( ABS(totfrac - un) > epsilocal) THEN
3762	WRITE(str1,'("This occurs on grid box", I8)') ji
3763	WRITE(str2,'("The sum over all fraction and error are ", E14.4, E14.4)') totfrac, ABS(totfrac - un)
3764	CALL ipslerr_p (3,'slowproc_checkveget', &
3765	"veget_max + frac_nobio is not equal to 1.", str1, str2)
3766	WRITE(,) "EPSILON =", epsilocal
3767	ENDIF
3768	ENDDO
3769
3770	!! 3.0 Verify that veget is smaller or equal to veget_max
3771	!!
3772	DO ji=1,nbpt
3773	DO jv=1,nvm
3774	IF ( jv == ibare_sechiba ) THEN
3775	IF ( ABS(veget(ji,jv) - veget_max(ji,jv)) > epsilocal ) THEN
3776	WRITE(str1,'("This occurs on grid box", I8)') ji
3777	WRITE(str2,'("The difference is ", E14.4)') veget(ji,jv) - veget_max(ji,jv)
3778	CALL ipslerr_p (3,'slowproc_checkveget', &
3779	"veget is not equal to veget_max on bare soil.", str1, str2)
3780	ENDIF
3781	ELSE
3782	IF ( veget(ji,jv) > veget_max(ji,jv) ) THEN
3783	WRITE(str1,'("This occurs on grid box", I8)') ji
3784	WRITE(str2,'("The values for veget and veget_max :", F8.4, F8.4)') veget(ji,jv), veget_max(ji,jv)
3785	CALL ipslerr_p (3,'slowproc_checkveget', &
3786	"veget is greater than veget_max.", str1, str2)
3787	ENDIF
3788	ENDIF
3789	ENDDO
3790	ENDDO
3791
3792	!! 4.0 Test tot_bare_soil in relation to the other variables
3793	!!
3794	DO ji=1,nbpt
3795	totfrac = zero
3796	DO jv=1,nvm
3797	totfrac = totfrac + (veget_max(ji,jv) - veget(ji,jv))
3798	ENDDO
3799	! add the bare soil fraction to totfrac
3800	totfrac = totfrac + veget(ji,ibare_sechiba)
3801	! do the test
3802	IF ( ABS(totfrac - tot_bare_soil(ji)) > epsilocal ) THEN
3803	WRITE(str1,'("This occurs on grid box", I8)') ji
3804	WRITE(str2,'("The values for tot_bare_soil, tot frac and error :", F8.4, F8.4, E14.4)') &
3805	& tot_bare_soil(ji), totfrac, ABS(totfrac - tot_bare_soil(ji))
3806	CALL ipslerr_p (3,'slowproc_checkveget', &
3807	"tot_bare_soil does not correspond to the total bare soil fraction.", str1, str2)
3808	ENDIF
3809	ENDDO
3810
3811	!! 5.0 Test that soiltile has the right sum
3812	!!
3813	DO ji=1,nbpt
3814	totfrac = SUM(soiltile(ji,:))
3815	IF ( ABS(totfrac - un) > epsilocal ) THEN
3816	WRITE(numout,*) "soiltile does not sum-up to one. This occurs on grid box", ji
3817	WRITE(numout,*) "The soiltile for ji are :", soiltile(ji,:)
3818	CALL ipslerr_p (2,'slowproc_checkveget', &
3819	"soiltile does not sum-up to one.", "", "")
3820	ENDIF
3821	ENDDO
3822
3823	END SUBROUTINE slowproc_checkveget
3824
3825
3826	!! ================================================================================================================================
3827	!! SUBROUTINE : slowproc_change_frac
3828	!!
3829	!>\BRIEF Update the vegetation fractions
3830	!!
3831	!! DESCRIPTION : Update the vegetation fractions. This subroutine is called in the same time step as lcchange in stomatelpj has
3832	!! has been done. This subroutine is called after the diagnostics have been written in sechiba_main.
3833	!!
3834	!! RECENT CHANGE(S): None
3835	!!
3836	!! MAIN OUTPUT VARIABLE(S): :: veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile
3837	!!
3838	!! REFERENCE(S) : None
3839	!!
3840	!! FLOWCHART : None
3841	!! \n
3842	!_ ================================================================================================================================
3843
3844	SUBROUTINE slowproc_change_frac(kjpindex, lai, &
3845	veget_max, veget, frac_nobio, totfrac_nobio, tot_bare_soil, soiltile, fraclut, nwdFraclut)
3846	!
3847	! 0. Declarations
3848	!
3849	! 0.1 Input variables
3850	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size - terrestrial pixels only
3851	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(in) :: lai !! Leaf area index (m^2 m^{-2})
3852
3853	! 0.2 Output variables
3854	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget_max !! Maximum fraction of vegetation type in the mesh (unitless)
3855	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(out) :: veget !! Fraction of vegetation type in the mesh (unitless)
3856	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT(out) :: frac_nobio !! Fraction of ice, lakes, cities etc. in the mesh
3857	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in the mesh
3858	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_bare_soil !! Total evaporating bare soil fraction in the mesh
3859	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out) :: soiltile !! Fraction of each soil tile within vegtot (0-1, unitless)
3860	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: fraclut !! Fraction of each landuse tile (0-1, unitless)
3861	REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out) :: nwdfraclut !! Fraction of non woody vegetation in each landuse tile (0-1, unitless)
3862
3863	! 0.3 Local variables
3864	INTEGER(i_std) :: ji, jv !! Loop index
3865
3866
3867	!! Update vegetation fractions with the values coming from the vegetation file read in slowproc_readvegetmax.
3868	!! Partial update has been taken into account for the case with DGVM and AGRICULTURE in slowproc_readvegetmax.
3869	veget_max = veget_max_new
3870	frac_nobio = frac_nobio_new
3871
3872	!! Verification and correction on veget_max, calculation of veget and soiltile.
3873	CALL slowproc_veget (kjpindex, lai, frac_nobio, totfrac_nobio, veget_max, veget, soiltile, fraclut, nwdFraclut)
3874
3875	!! Calculate tot_bare_soil needed in hydrol, diffuco and condveg (fraction of bare soil in the mesh)
3876	tot_bare_soil(:) = veget_max(:,1)
3877	DO jv = 2, nvm
3878	DO ji =1, kjpindex
3879	tot_bare_soil(ji) = tot_bare_soil(ji) + (veget_max(ji,jv) - veget(ji,jv))
3880	ENDDO
3881	END DO
3882
3883	!! Do some basic tests on the surface fractions updated above
3884	CALL slowproc_checkveget(kjpindex, frac_nobio, veget_max, veget, tot_bare_soil, soiltile)
3885
3886	END SUBROUTINE slowproc_change_frac
3887
3888	END MODULE slowproc

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

Download in other formats: