Context Navigation

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

Last change on this file since 7325 was 7325, checked in by josefine.ghattas, 3 years ago
Add test for each of the variables njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, clayfraction, sandfraction if they were not found in restart file. If at least one was missing, activate call to slowproc_soilt. This commit is needed if using old restart file were clayfraction was in the restart file but not the 7 new variables (ks, nvan, avan, mcr, mcs, mcfc, mcw) which were added with modifications from branch SP-MIP. See ticket #800
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 208.7 KB

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

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

Download in other formats: