Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

slowproc.f90 @ 7541

Last change on this file since 7541 was 7541, checked in by fabienne.maignan, 2 years ago
Zhang publication on coupling factor
File size: 207.5 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_2.2_r7266/ORCHIDEE/src_sechiba/slowproc.f90 @ 7541

Download in other formats: