Context Navigation

source: branches/publications/ORCHIDEE_CN_CAN_r5698/src_sechiba/slowproc.f90 @ 7346

Last change on this file since 7346 was 5691, checked in by sebastiaan.luyssaert, 6 years ago
DEV: tested with 13, 37 and 64 PFTs with LCC on different pixels. Some configuration run for 20 years on a given pixel, other crash on another pixel. There is a mass balance problem in sapiens_lcc (ticket #482). This commit fixes a problem with PFT1 in littercalc. This PFT is now fully integrated in LCC and subsequent litter and soil dynamics. veget_max was changed in veget_cov_max where appropriate, a typo in enerbil was corrected.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 243.8 KB

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

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

Download in other formats: