Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

slowproc.f90 @ 7329

Last change on this file since 7329 was 7191, checked in by haicheng.zhang, 3 years ago
Implementing latral transfers of sediment and POC from land to ocean through river in ORCHILEAK
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 240.3 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-Clateral/src_sechiba/slowproc.f90 @ 7329

Download in other formats: