Context Navigation

slowproc.f90 @ 6319

Last change on this file since 6319 was 6319, checked in by josefine.ghattas, 5 years ago

Improvment for the ESM CO2 configuration:

Separate variable fco2_lu into 3 parts: fco2_lu, fco2_wh and fco2_ha
Move calculation of co2_flux from dt_sechiba time-step to daily time-step (in the part for do_slow)
Removed co2_flux and fco2_lu from stomate_intialize argument list. These variables were never used in the intialization phase.
Add co2_flux, and fco2_wh, fco2_ha to restart file
Corrected output unit for nee to be consistent with LMDZ and stomate output variables. It is now in kgC/m2/s.
Corrected output for znetco2
Added fCO2_fWoodharvest and fCO2_fHarvest as new possible tracers in LMDZ (intersurf).
Added diagnostic output for fCO2_fWoodharvest and fCO2_fHarvest

Cadule

Property svn:keywords set to HeadURL Date Author Revision

File size: 184.9 KB

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

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

Context Navigation

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

Download in other formats: