Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_highres.f90 @ 8504

Last change on this file since 8504 was 8504, checked in by josefine.ghattas, 7 weeks ago
Integrated correction done in the trunk [8503] for compilation without XIOS.
Property svn:executable set to ``*
File size: 243.5 KB

Line
1	! ================================================================================================================================
2	! MODULE : routing_highres
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module routes the water over the continents into the oceans and computes the water
10	!! stored in floodplains or taken for irrigation.
11	!!
12	!!\n DESCRIPTION: None
13	!!
14	!! RECENT CHANGE(S): Now works together with the routing pre-processor : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp
15	!!
16	!! REFERENCE(S) :
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-ROUTING/ORCHIDEE/src_sechiba/routing.f90 $
20	!! $Date: 2022-03-24 11:25:05 +0100 (Do, 24 MÃ€r 2022) $
21	!! $Revision: 7545 $
22	!! \n
23	!_ ================================================================================================================================
24	!
25	!
26	! Histoire Salee
27	!---------------
28	! La douce riviere
29	! Sortant de son lit
30	! S'est jetee ma chere
31	! dans les bras mais oui
32	! du beau fleuve
33	!
34	! L'eau coule sous les ponts
35	! Et puis les flots s'emeuvent
36	! - N'etes vous pas au courant ?
37	! Il parait que la riviere
38	! Va devenir mer
39	! Roland Bacri
40	!
41
42
43	MODULE routing_highres
44
45	USE ioipsl
46	USE xios_orchidee
47	USE ioipsl_para
48	USE constantes
49	USE constantes_var
50	USE constantes_soil
51	USE pft_parameters
52	USE sechiba_io_p
53	USE interpol_help
54	USE grid
55	USE mod_orchidee_para
56
57	USE haversine
58
59	IMPLICIT NONE
60	PRIVATE
61	PUBLIC :: routing_highres_main, routing_highres_initialize, routing_highres_finalize, routing_highres_clear, routing_highres_xios_initialize
62
63	INTERFACE routing_hr_landgather
64	MODULE PROCEDURE routing_hr_landgather_i1, routing_hr_landgather_i2, routing_hr_landgather_r
65	END INTERFACE routing_hr_landgather
66
67	INTEGER(i_std),PARAMETER :: WaterCp=1000.*4.1813 !! water heat capacity in J/Kg/K
68
69	!! PARAMETERS
70	INTEGER(i_std), SAVE :: nbasmax=-1 !! The maximum number of basins we wish to have per grid box (truncation of the model) (unitless)
71	INTEGER(i_std), SAVE :: nbasmon = 4 !! Number of basins to be monitored
72	INTEGER(i_std), SAVE :: inflows=-1 !! The maximum number of inflows (unitless)
73	INTEGER(i_std), SAVE :: nbvmax !! The maximum number of basins we can handle at any time during the generation of the maps (unitless)
74	!$OMP THREADPRIVATE(nbvmax)
75	REAL(r_std), SAVE :: fast_tcst = -1. !! Property of the fast reservoir, (s/km)
76	!$OMP THREADPRIVATE(fast_tcst)
77	REAL(r_std), SAVE :: slow_tcst = -1. !! Property of the slow reservoir, (s/km)
78	!$OMP THREADPRIVATE(slow_tcst)
79	REAL(r_std), SAVE :: stream_tcst = -1. !! Property of the stream reservoir, (s/km)
80	!$OMP THREADPRIVATE(stream_tcst)
81	REAL(r_std), SAVE :: flood_tcst = -1. !! Property of the floodplains reservoir, (s/km)
82	!$OMP THREADPRIVATE(flood_tcst)
83	REAL(r_std), SAVE :: swamp_cst = -1. !! Fraction of the river transport that flows to the swamps (unitless;0-1)
84	!$OMP THREADPRIVATE(swamp_cst)
85	REAL(r_std), SAVE :: lim_floodcri = -1. !! Minimal orog diff between two consecutive floodplains htu (m)
86	!$OMP THREADPRIVATE(lim_floodcri)
87	!
88	! Relation between volume and fraction of floodplains
89	!
90	REAL(r_std), SAVE :: betap = 0.5 !! Ratio of the basin surface intercepted by ponds and the maximum surface of ponds (unitless;0-1)
91	!$OMP THREADPRIVATE(betap)
92	REAL(r_std), SAVE :: rfloodmax = 0.5 !! Maximal discharge reducer when there are floodplains
93	!$OMP THREADPRIVATE(rfloodmax)
94	REAL(r_std), SAVE :: overflow_tcst = 5 !! Maximal discharge reducer when there are floodplains
95	!$OMP THREADPRIVATE(overflow_tcst)
96	INTEGER(i_std), SAVE :: overflow_repetition = 1 !! Number of repetition of overflow for each routing step
97	!$OMP THREADPRIVATE(overflow_repetition)
98	! Soil temperature depth to be used to estimate runoff and drainage temperatures
99	!
100	REAL(r_std), PARAMETER, DIMENSION(2) :: runofftempdepth = (/ 0.0, 0.3 /) !! Layer which will determine the temperature of runoff
101	REAL(r_std), PARAMETER, DIMENSION(2) :: drainagetempdepth = (/ 3.0, 90.0 /) !! Layer which will determine the temperature of runoff
102	!
103	!
104	! Relation between maximum surface of ponds and basin surface, and drainage (mm/j) to the slow_res
105	!
106	REAL(r_std), PARAMETER :: pond_bas = 50.0 !! [DISPENSABLE] - not used
107	REAL(r_std), SAVE :: pondcri = 2000.0 !! Potential height for which all the basin is a pond (mm)
108	!$OMP THREADPRIVATE(pondcri)
109	!
110	REAL(r_std), PARAMETER :: maxevap_lake = 7.5/86400. !! Maximum evaporation rate from lakes (kg/m^2/s)
111	!
112	REAL(r_std),SAVE :: dt_routing !! Routing time step (s)
113	!$OMP THREADPRIVATE(dt_routing)
114	!
115	INTEGER(i_std), SAVE :: ntemp_layer = 4 !! Number of layers to be taken to determine the ground water temperature.
116	!$OMP THREADPRIVATE(ntemp_layer)
117	INTEGER(i_std), SAVE :: diagunit = 87 !! Diagnostic file unit (unitless)
118	!$OMP THREADPRIVATE(diagunit)
119	!
120	! Logicals to control model configuration
121	!
122	LOGICAL, SAVE :: dofloodinfilt = .FALSE. !! Logical to choose if floodplains infiltration is activated or not (true/false)
123	!$OMP THREADPRIVATE(dofloodinfilt)
124	LOGICAL, SAVE :: dofloodoverflow = .FALSE. !! Logical to choose if floodplains overflow is activated or not (true/false)
125	!$OMP THREADPRIVATE(dofloodoverflow)
126	LOGICAL, SAVE :: doswamps = .FALSE. !! Logical to choose if swamps are activated or not (true/false)
127	!$OMP THREADPRIVATE(doswamps)
128	LOGICAL, SAVE :: doponds = .FALSE. !! Logical to choose if ponds are activated or not (true/false)
129	!$OMP THREADPRIVATE(doponds)
130	REAL(r_std), SAVE :: conduct_factor = 1. !! Adjustment factor for floodplains reinfiltration
131	!$OMP THREADPRIVATE(conduct_factor)
132	!
133	! The variables describing the basins and their routing, need to be in the restart file.
134	!
135	INTEGER(i_std), SAVE :: num_largest = 200 !! Number of largest river basins which should be treated as independently as rivers
136	!! (not flow into ocean as diffusion coastal flow) (unitless)
137	!$OMP THREADPRIVATE(num_largest)
138	!
139	CHARACTER(LEN=80),SAVE :: graphfilename="routing_graph.nc"
140	!$OMP THREADPRIVATE(graphfilename)
141	REAL(r_std), SAVE :: undef_graphfile
142	!$OMP THREADPRIVATE(undef_graphfile)
143	REAL(r_std), SAVE :: graphfile_version = 0.0
144	!$OMP THREADPRIVATE(graphfile_version)
145	REAL(r_std), SAVE :: maxtimestep = 1800.0 !! A reasonalble maximum time step. Actual value to be read from graphfile.
146	!$OMP THREADPRIVATE(maxtimestep)
147	REAL(r_std), SAVE :: time_counter !! Time counter (s)
148	!$OMP THREADPRIVATE(time_counter)
149	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_loc !! Surface of basin (m^2)
150	!$OMP THREADPRIVATE(routing_area_loc)
151	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_loc !! Topographic index of the retention time (m)
152	!$OMP THREADPRIVATE(topo_resid_loc)
153	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: stream_resid_loc !! Topographic index of the retention time (m)
154	!$OMP THREADPRIVATE(stream_resid_loc)
155	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_loc !! Grid into which the basin flows (unitless)
156	!$OMP THREADPRIVATE(route_togrid_loc)
157	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_loc !! Basin in to which the water goes (unitless)
158	!$OMP THREADPRIVATE(route_tobasin_loc)
159	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_loc !! Number of basin into current one (unitless)
160	!$OMP THREADPRIVATE(route_nbintobas_loc)
161	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_loc !! ID of basin (unitless)
162	!$OMP THREADPRIVATE(global_basinid_loc)
163	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydrodiag_loc !! Variable to diagnose the hydrographs
164	!$OMP THREADPRIVATE(hydrodiag_loc)
165	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: HTUdiag_loc !! Variable to diagnose the hydrographs
166	!$OMP THREADPRIVATE(HTUdiag_loc)
167	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: HTUdiag_glo !! Variable to diagnose the hydrographs
168	!$OMP THREADPRIVATE(HTUdiag_glo)
169	LOGICAL, SAVE :: MonitoringinGraph=.FALSE.
170	LOGICAL, SAVE :: ReadGraph=.FALSE.
171	LOGICAL, SAVE :: ReadMonitoring=.FALSE.
172	REAL(r_std), SAVE :: stream_maxresid
173	!$OMP THREADPRIVATE(stream_maxresid)
174	!
175	! parallelism
176	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
177	!$OMP THREADPRIVATE(routing_area_glo)
178	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m)
179	!$OMP THREADPRIVATE(topo_resid_glo)
180	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: stream_resid_glo !! Topographic index of the retention time (m)
181	!$OMP THREADPRIVATE(stream_resid_glo)
182	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless)
183	!$OMP THREADPRIVATE(route_togrid_glo)
184	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless)
185	!$OMP THREADPRIVATE(route_tobasin_glo)
186	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_glo !! Number of basin into current one (unitless)
187	!$OMP THREADPRIVATE(route_nbintobas_glo)
188	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless)
189	!$OMP THREADPRIVATE(global_basinid_glo)
190	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydrodiag_glo !! Variable to diagnose the hydrographs
191	!$OMP THREADPRIVATE(hydrodiag_glo)
192	!
193	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: routing_area !! Surface of basin (m^2)
194	!$OMP THREADPRIVATE(routing_area)
195	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: topo_resid !! Topographic index of the retention time (m)
196	!$OMP THREADPRIVATE(topo_resid)
197	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: stream_resid !! Topographic index of the retention time (m)
198	!$OMP THREADPRIVATE(stream_resid)
199	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_togrid !! Grid into which the basin flows (unitless)
200	!$OMP THREADPRIVATE(route_togrid)
201	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_tobasin !! Basin in to which the water goes (unitless)
202	!$OMP THREADPRIVATE(route_tobasin)
203	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_nbintobas !! Number of basin into current one (unitless)
204	!$OMP THREADPRIVATE(route_nbintobas)
205	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: global_basinid !! ID of basin (unitless)
206	!$OMP THREADPRIVATE(global_basinid)
207	INTEGER(i_std), SAVE, POINTER, DIMENSION(:) :: hydrodiag !! Variable to diagnose the hydrographs
208	!$OMP THREADPRIVATE(hydrodiag)
209	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slowflow_diag !! Diagnostic slow flow hydrographs (kg/dt)
210	!$OMP THREADPRIVATE(slowflow_diag)
211	!
212	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigated !! Area equipped for irrigation in each grid box (m^2)
213	!$OMP THREADPRIVATE(irrigated)
214	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodplains_glo !! Maximal surface which can be inundated in each grid box (m^2)
215	!$OMP THREADPRIVATE(floodplains_glo)
216	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodplains_loc !! Maximal surface which can be inundated in each grid box (m^2)
217	!$OMP THREADPRIVATE(floodplains_loc)
218	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: floodplains !! Maximal surface which can be inundated in each grid box (m^2)
219	!$OMP THREADPRIVATE(floodplains)
220	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodmap
221	!! Floodplains Fraction for each grid point.
222	!$OMP THREADPRIVATE(floodmap)
223
224	!!!
225
226	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: tempdiag_mean !! Averaged soil temperatures
227	!$OMP THREADPRIVATE(tempdiag_mean)
228	!
229	! FLOOD OVERFLOW
230	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: orog_min_glo !!
231	!$OMP THREADPRIVATE(orog_min_glo)!
232	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: orog_min_loc !!
233	!$OMP THREADPRIVATE(orog_min_loc)
234	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: orog_min !!
235	!$OMP THREADPRIVATE(orog_min)
236	!
237	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_innum_glo !!
238	!$OMP THREADPRIVATE(route_innum_glo)
239	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_innum_loc !!
240	!$OMP THREADPRIVATE(route_innum_loc)
241	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_innum !!
242	!$OMP THREADPRIVATE(route_innum)
243	!
244	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_ingrid_glo !!
245	!$OMP THREADPRIVATE(route_ingrid_glo)
246	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_ingrid_loc !!
247	!$OMP THREADPRIVATE(route_ingrid_loc)
248	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:,:) :: route_ingrid !!
249	!$OMP THREADPRIVATE(route_ingrid)
250	!
251	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_inbasin_glo !!
252	!$OMP THREADPRIVATE(route_inbasin_glo)
253	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_inbasin_loc !!
254	!$OMP THREADPRIVATE(route_inbasin_loc)
255	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:,:) :: route_inbasin !!
256	!$OMP THREADPRIVATE(route_inbasin)
257	!
258	!!!
259	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: swamp !! Maximal surface of swamps in each grid box (m^2)
260	!$OMP THREADPRIVATE(swamp)
261	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: fp_beta_glo !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless)
262	!$OMP THREADPRIVATE(fp_beta_glo)
263	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: fp_beta_loc !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless)
264	!$OMP THREADPRIVATE(fp_beta_loc)
265	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: fp_beta !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless)
266	!$OMP THREADPRIVATE(fp_beta)
267	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodcri_glo !! Potential height for which all the basin is a pond (mm)
268	!$OMP THREADPRIVATE(floodcri_glo)
269	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodcri_loc !! Potential height for which all the basin is a pond (mm)
270	!$OMP THREADPRIVATE(floodcri_loc)
271	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: floodcri !! Potential height for which all the basin is a pond (mm)
272	!$OMP THREADPRIVATE(floodcri)
273	!
274	! The reservoirs, also to be put into the restart file.
275	!
276	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_reservoir !! Water amount in the fast reservoir (kg)
277	!$OMP THREADPRIVATE(fast_reservoir)
278	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: slow_reservoir !! Water amount in the slow reservoir (kg)
279	!$OMP THREADPRIVATE(slow_reservoir)
280	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: stream_reservoir !! Water amount in the stream reservoir (kg)
281	!$OMP THREADPRIVATE(stream_reservoir)
282	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir !! Water amount in the floodplains reservoir (kg)
283	!$OMP THREADPRIVATE(flood_reservoir)
284	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_reservoir !! Water amount in the lake reservoir (kg)
285	!$OMP THREADPRIVATE(lake_reservoir)
286	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_reservoir !! Water amount in the pond reservoir (kg)
287	!$OMP THREADPRIVATE(pond_reservoir)
288	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas !! Flooded fraction per basin (unitless;0-1)
289	!$OMP THREADPRIVATE(flood_frac_bas)
290	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_frac !! Pond fraction per grid box (unitless;0-1)
291	!$OMP THREADPRIVATE(pond_frac)
292	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_height !! Floodplain height (mm)
293	!$OMP THREADPRIVATE(flood_height)
294	!
295	! Reservoir temperatures
296	!
297	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_temp !! Water temperature in the fast reservoir (K)
298	!$OMP THREADPRIVATE(fast_temp)
299	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: slow_temp !! Water temperature in the slow reservoir (K)
300	!$OMP THREADPRIVATE(slow_temp)
301	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: stream_temp !! Water temperature in the stream reservoir (K)
302	!$OMP THREADPRIVATE(stream_temp)
303	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: streamlimit !!
304	!$OMP THREADPRIVATE(streamlimit)
305	!
306	! The accumulated fluxes.
307	!
308	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodout_mean !! Accumulated flow out of floodplains (kg/m^2/dt)
309	!$OMP THREADPRIVATE(floodout_mean)
310	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: runoff_mean !! Accumulated runoff (kg/m^2/dt)
311	!$OMP THREADPRIVATE(runoff_mean)
312	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: drainage_mean !! Accumulated drainage (kg/m^2/dt)
313	!$OMP THREADPRIVATE(drainage_mean)
314	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: transpot_mean !! Mean potential transpiration from the plants (kg/m^2/dt)
315	!$OMP THREADPRIVATE(transpot_mean)
316	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: precip_mean !! Accumulated precipitation (kg/m^2/dt)
317	!$OMP THREADPRIVATE(precip_mean)
318	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: humrel_mean !! Mean soil moisture stress, mean root extraction potential (unitless)
319	!$OMP THREADPRIVATE(humrel_mean)
320	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: totnobio_mean !! Mean last total fraction of no bio (unitless;0-1)
321	!$OMP THREADPRIVATE(totnobio_mean)
322	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: vegtot_mean !! Mean potentially vegetated fraction (unitless;0-1)
323	!$OMP THREADPRIVATE(vegtot_mean)
324	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: k_litt_mean !! Mean averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
325	!$OMP THREADPRIVATE(k_litt_mean)
326	!
327	! The averaged outflow fluxes.
328	!
329	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lakeinflow_mean !! Mean lake inflow (kg/m^2/dt)
330	!$OMP THREADPRIVATE(lakeinflow_mean)
331	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: returnflow_mean !! Mean water flow from lakes and swamps which returns to the grid box.
332	!! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt)
333	!$OMP THREADPRIVATE(returnflow_mean)
334	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: reinfiltration_mean !! Mean water flow which returns to the grid box (kg/m^2/dt)
335	!$OMP THREADPRIVATE(reinfiltration_mean)
336	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigation_mean !! Mean irrigation flux.
337	!! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
338	!$OMP THREADPRIVATE(irrigation_mean)
339	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean !! Mean Outflow of the major rivers.
340	!! The flux will be located on the continental grid but this should be a coastal point (kg/dt)
341	!$OMP THREADPRIVATE(riverflow_mean)
342	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: coastalflow_mean !! Mean outflow on coastal points by small basins.
343	!! This is the water which flows in a disperse way into the ocean (kg/dt)
344	!$OMP THREADPRIVATE(coastalflow_mean)
345	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodtemp !! Temperature to decide if floodplains work (K)
346	!$OMP THREADPRIVATE(floodtemp)
347	INTEGER(i_std), SAVE :: floodtemp_lev !! Temperature level to decide if floodplains work (K)
348	!$OMP THREADPRIVATE(floodtemp_lev)
349	!
350	! Diagnostic variables ... well sort of !
351	!
352	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_netereq !! Irrigation requirement (water requirements by the crop for its optimal growth (kg/m^2/dt)
353	!$OMP THREADPRIVATE(irrig_netereq)
354	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrographs !! Hydrographs at the outflow of the grid box for major basins (kg/dt)
355	!$OMP THREADPRIVATE(hydrographs)
356	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrotemp !! Temperature of the largest river (in the HTUdiag sense) in the grid (K)
357	!$OMP THREADPRIVATE(hydrotemp)
358	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUhgmon !! Hydrographs to be monitored on specific HTUs (kg/dt)
359	!$OMP THREADPRIVATE(HTUhgmon)
360	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUhgmon_glo !! Hydrographs to be monitored on specific HTUs (kg/dt)
361	!$OMP THREADPRIVATE(HTUhgmon_glo)
362	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUtempmon !! Temperature to be monitored on specific HTUs (K)
363	!$OMP THREADPRIVATE(HTUtempmon)
364	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUtempmon_glo !! Temperature to be monitored on specific HTUs (K)
365	!$OMP THREADPRIVATE(HTUtempmon_glo)
366	!
367	! Diagnostics for the various reservoirs we use (Kg/m^2)
368	!
369	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: fast_diag !! Diagnostic for the fast reservoir (kg/m^2)
370	!$OMP THREADPRIVATE(fast_diag)
371	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slow_diag !! Diagnostic for the slow reservoir (kg/m^2)
372	!$OMP THREADPRIVATE(slow_diag)
373	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: stream_diag !! Diagnostic for the stream reservoir (kg/m^2)
374	!$OMP THREADPRIVATE(stream_diag)
375	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_diag !! Diagnostic for the floodplain reservoir (kg/m^2)
376	!$OMP THREADPRIVATE(flood_diag)
377	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_diag !! Diagnostic for the pond reservoir (kg/m^2)
378	!$OMP THREADPRIVATE(pond_diag)
379	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_diag !! Diagnostic for the lake reservoir (kg/m^2)
380	!$OMP THREADPRIVATE(lake_diag)
381
382	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: mask_coast !! Mask with coastal gridcells on local grid(1/0)
383	!$OMP THREADPRIVATE(mask_coast)
384	REAL(r_std), SAVE :: max_lake_reservoir !! Maximum limit of water in lake_reservoir [kg/m2]
385	!$OMP THREADPRIVATE(max_lake_reservoir)
386	INTEGER(i_std), SAVE :: nb_coast_gridcells !! Number of gridcells which can receive coastalflow
387	!$OMP THREADPRIVATE(nb_coast_gridcells)
388
389
390	CONTAINS
391	!! =============================================================================================================================
392	!! SUBROUTINE: routing_highres_initialize
393	!!
394	!>\BRIEF Initialize the routing module
395	!!
396	!! DESCRIPTION: Initialize the routing module. Read from restart file or read the routing.nc file to initialize the
397	!! routing scheme.
398	!!
399	!! RECENT CHANGE(S)
400	!!
401	!! REFERENCE(S)
402	!!
403	!! FLOWCHART
404	!! \n
405	!_ ==============================================================================================================================
406
407	SUBROUTINE routing_highres_initialize( kjit, nbpt, index, &
408	rest_id, hist_id, hist2_id, lalo, &
409	neighbours, resolution, contfrac, tempdiag, &
410	returnflow, reinfiltration, irrigation, riverflow, &
411	coastalflow, flood_frac, flood_res )
412
413	IMPLICIT NONE
414
415	!! 0.1 Input variables
416	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
417	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
418	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless)
419	INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless)
420	INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless)
421	INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless)
422	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
423
424	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb) !! Vector of neighbours for each grid point
425	!! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
426	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m)
427	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
428	REAL(r_std), INTENT(in) :: tempdiag(nbpt,ngrnd) !! Diagnostic soil temperature profile
429
430	!! 0.2 Output variables
431	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box.
432	!! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt)
433	REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
434	REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
435	REAL(r_std), INTENT(out) :: riverflow(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
436
437	REAL(r_std), INTENT(out) :: coastalflow(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
438	REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
439	REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
440
441	!! 0.3 Local variables
442	REAL(r_std), DIMENSION(nbp_glo) :: mask_coast_glo !! Mask with coastal gridcells on global grid (1/0)
443	LOGICAL :: init_irrig !! Logical to initialize the irrigation (true/false)
444	LOGICAL :: init_flood !! Logical to initialize the floodplains (true/false)
445	LOGICAL :: init_swamp !! Logical to initialize the swamps (true/false)
446	INTEGER :: ig, ib, rtg, rtb !! Index
447	REAL(r_std) :: stream_tcst_orig
448	INTEGER :: ier !! Error handeling
449	!_ ================================================================================================================================
450
451	!
452	! do initialisation
453	!
454	nbvmax = 440
455	! Here we will allocate the memory and get the fixed fields from the restart file.
456	! If the info is not found then we will compute the routing map.
457	!
458
459	CALL routing_hr_init (kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
460	riverflow, coastalflow, flood_frac, flood_res, tempdiag, rest_id)
461
462	routing_area => routing_area_loc
463	floodplains => floodplains_loc
464	topo_resid => topo_resid_loc
465	stream_resid => stream_resid_loc
466	route_togrid => route_togrid_loc
467	route_tobasin => route_tobasin_loc
468	global_basinid => global_basinid_loc
469	hydrodiag => hydrodiag_loc
470	fp_beta => fp_beta_loc
471	floodcri => floodcri_loc
472	!
473	route_innum => route_innum_loc
474	route_ingrid => route_ingrid_loc
475	route_inbasin => route_inbasin_loc
476	orog_min => orog_min_loc
477
478	! This routine computes the routing map if the route_togrid_glo is undefined. This means that the
479	! map has not been initialized during the restart process..
480	!
481	!! Reads in the map of the basins and flow directions to construct the catchments of each grid box
482	!
483	IF ( ReadGraph .OR. ReadMonitoring) THEN
484	CALL routing_hr_basins_p(nbpt, lalo, neighbours, resolution, contfrac)
485	ENDIF
486	! Keep the information so we can check the time step.
487	stream_tcst_orig = stream_tcst
488	!
489	IF (stream_tcst .LE. 0 .OR. fast_tcst .LE. 0 .OR. slow_tcst .LE. 0 .OR. flood_tcst .LE. 0 ) THEN
490	CALL ipslerr(3,'routing_highres_initialize',' The time constants of the routing reservoirs were not initialized. ', &
491	'Please check if they are present in the HTU graph file', ' ')
492	ELSE
493	!
494	!> The time constants for the various reservoirs should be read from the graph file
495	!> produced by routingpp (https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp). They are
496	!> also saved in the restart file so that we do not need to read the graph file at each restart.
497	!> But once they are set in the model the user can changed them through the run.def.
498	!> This is a useful option to test values but should not be an operational solution. The
499	!> correct value should be given to the model through the graph file.
500	!> The getin_p operation cannot be done earlier as in routing_hr_init above these constant
501	!> might not yet be known.
502	!
503	!Config Key = SLOW_TCST
504	!Config Desc = Time constant for the slow reservoir
505	!Config If = RIVER_ROUTING
506	!Config Def = 25.0
507	!Config Help = This parameters allows the user to fix the
508	!Config time constant (s/km) of the slow reservoir
509	!Config in order to get better river flows for
510	!Config particular regions.
511	!Config Units = [days]
512	!
513	CALL getin_p('SLOW_TCST', slow_tcst)
514	!
515	!Config Key = FAST_TCST
516	!Config Desc = Time constant for the fast reservoir
517	!Config If = RIVER_ROUTING
518	!Config Def = 3.0
519	!Config Help = This parameters allows the user to fix the
520	!Config time constant (s/km) of the fast reservoir
521	!Config in order to get better river flows for
522	!Config particular regions.
523	!Config Units = [days]
524	CALL getin_p('FAST_TCST', fast_tcst)
525
526	!Config Key = STREAM_TCST
527	!Config Desc = Time constant for the stream reservoir
528	!Config If = RIVER_ROUTING
529	!Config Def = 0.24
530	!Config Help = This parameters allows the user to fix the
531	!Config time constant (s/km) of the stream reservoir
532	!Config in order to get better river flows for
533	!Config particular regions.
534	!Config Units = [days]
535	CALL getin_p('STREAM_TCST', stream_tcst)
536
537	!Config Key = FLOOD_TCST
538	!Config Desc = Time constant for the flood reservoir
539	!Config If = RIVER_ROUTING
540	!Config Def = 4.0
541	!Config Help = This parameters allows the user to fix the
542	!Config time constant (s/km) of the flood reservoir
543	!Config in order to get better river flows for
544	!Config particular regions.
545	!Config Units = [days]
546	CALL getin_p('FLOOD_TCST', flood_tcst)
547
548	!Config Key = SWAMP_CST
549	!Config Desc = Fraction of the river that flows back to swamps
550	!Config If = RIVER_ROUTING
551	!Config Def = 0.2
552	!Config Help = This parameters allows the user to fix the
553	!Config fraction of the river transport
554	!Config that flows to swamps
555	!Config Units = [-]
556	CALL getin_p('SWAMP_CST', swamp_cst)
557	!
558	!Config Key = LIM_FLOODCRI
559	!Config Desc = Difference of orography between floodplains HTUs.
560	!Config If = RIVER_ROUTING
561	!Config Def = 0.3
562	!Config Help = This parameters allows the user to fix the
563	!Config minimal difference of orography between two consecutive
564	!Config floodplains HTU.
565	!Config Units = [meter]
566	CALL getin_p('LIM_FLOODCRI', lim_floodcri)
567	!
568	ENDIF
569	!
570	! Verify that the time step is compatible with the graph file.
571	! If the user has changed the time constant of the stream reservoir then
572	! the maximum time step needs to be adjusted.
573	!
574	IF ( stream_tcst_orig == 0 ) THEN
575	WRITE(,) "routing_highres_initialize : Update stream_tcst ", stream_tcst_orig, stream_tcst
576	stream_tcst_orig = stream_tcst
577	ENDIF
578	IF ( dt_routing > maxtimestep/stream_tcst_orig*stream_tcst ) THEN
579	WRITE(,) "routing_highres_initialize : Chosen time step : ", dt_routing
580	WRITE(,) "routing_highres_initialize : Recommended time step : ", maxtimestep/stream_tcst_orig*stream_tcst
581	CALL ipslerr_p(2,'routing_highres_initialize','The chosen time step is larger than the value recommended','in the graph file.','')
582	ENDIF
583	!
584	!
585	!
586	IF (dofloodoverflow) THEN
587	CALL routing_hr_inflows(nbp_glo, nbasmax, inflows, floodplains_glo,route_innum_glo,route_ingrid_glo,route_inbasin_glo)
588	END IF
589
590	!! Create a mask containing all possible coastal gridcells and count total number of coastal gridcells
591	IF (is_root_prc) THEN
592	mask_coast_glo(:)=0
593	DO ib=1,nbasmax
594	DO ig=1,nbp_glo
595	rtg = route_togrid_glo(ig,ib)
596	rtb = route_tobasin_glo(ig,ib)
597	! Coastal gridcells are stored in nbasmax+2
598	IF ( rtb == nbasmax+2) THEN
599	mask_coast_glo(rtg) = 1
600	END IF
601	END DO
602	END DO
603	nb_coast_gridcells=SUM(mask_coast_glo)
604	IF (printlev>=3) WRITE(numout,*) 'Number of coastal gridcells = ', nb_coast_gridcells
605
606	IF (nb_coast_gridcells == 0)THEN
607	CALL ipslerr(3,'routing_highres_initialize',&
608	'Number of coastal gridcells is zero for routing. ', &
609	'If this is a global run, this is an error.',&
610	'If this is a regional run, please check to make sure your region includes a full basin or turn routing off.')
611	ENDIF
612
613	ENDIF
614	CALL bcast(nb_coast_gridcells)
615
616	ALLOCATE(mask_coast(nbpt), stat=ier)
617	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_inititalize','Pb in allocate for mask_coast','','')
618	CALL scatter(mask_coast_glo, mask_coast)
619	!
620	! Do we have what we need if we want to do irrigation
621	!! Initialisation of flags for irrigated land, flood plains and swamps
622	!
623	init_irrig = .FALSE.
624	IF ( do_irrigation ) THEN
625	IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE.
626	END IF
627
628	init_flood = .FALSE.
629	IF ( do_floodplains ) THEN
630	IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) init_flood = .TRUE.
631	END IF
632
633	init_swamp = .FALSE.
634	IF ( doswamps ) THEN
635	IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE.
636	END IF
637
638	!! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree
639	!! base data set to the resolution of the model.
640
641	!IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN
642	! CALL routing_hr_irrigmap(nbpt, index, lalo, neighbours, resolution, &
643	! contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)
644	!ENDIF
645
646	IF (printlev >= 5) WRITE(numout,*) 'End of routing_highres_initialize'
647
648	END SUBROUTINE routing_highres_initialize
649
650	!! =============================================================================================================================
651	!! SUBROUTINE: routing_highres_xios_initialize
652	!!
653	!>\BRIEF Initialize xios dependant defintion before closing context defintion
654	!!
655	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion.
656	!! This subroutine is called before the xios context is closed.
657	!!
658	!! RECENT CHANGE(S): None
659	!!
660	!! REFERENCE(S): None
661	!!
662	!! FLOWCHART: None
663	!! \n
664	!_ ==============================================================================================================================
665
666	SUBROUTINE routing_highres_xios_initialize
667	USE xios_orchidee
668	IMPLICIT NONE
669
670	INTEGER(i_std) ::ib
671
672	!
673	! If the routing_graph file is available we will extract the information in the dimensions
674	! and parameters.
675	!
676	!Config Key = ROUTING_FILE
677	!Config Desc = Name of file which contains the routing information graph on the model grid
678	!Config If = RIVER_ROUTING
679	!Config Def = routing.nc
680	!Config Help = The file provided here should allows to route the water from one HTU
681	!Config to another. The RoutingPP code needs to be used in order to generate
682	!Config the routing graph for the model grid.
683	!Config More details on : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp
684	!Config Units = [FILE]
685	!
686	graphfilename = 'routing_graph.nc'
687	CALL getin('ROUTING_FILE',graphfilename)
688	CALL routing_hr_graphinfo(graphfilename, nbasmax, inflows, nbasmon, undef_graphfile, stream_tcst, fast_tcst, slow_tcst, &
689	& flood_tcst, swamp_cst, lim_floodcri)
690
691	CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/))
692	CALL xios_orchidee_addaxis("nbasmon", nbasmon, (/(REAL(ib,r_std),ib=1,nbasmon)/))
693
694	END SUBROUTINE routing_highres_xios_initialize
695
696	!! ================================================================================================================================
697	!! SUBROUTINE : routing_highres_main
698	!!
699	!>\BRIEF This module routes the water over the continents (runoff and
700	!! drainage produced by the hydrol module) into the oceans.
701	!!
702	!! DESCRIPTION (definitions, functional, design, flags):
703	!! The routing scheme (Polcher, 2003) carries the water from the runoff and drainage simulated by SECHIBA
704	!! to the ocean through reservoirs, with some delay. The routing scheme is based on
705	!! a parametrization of the water flow on a global scale (Miller et al., 1994; Hagemann
706	!! and Dumenil, 1998). Given the global map of the main watersheds (Oki et al., 1999;
707	!! Fekete et al., 1999; Vorosmarty et al., 2000) which delineates the boundaries of subbasins
708	!! and gives the eight possible directions of water flow within the pixel, the surface
709	!! runoff and the deep drainage are routed to the ocean. The time-step of the routing is one day.
710	!! The scheme also diagnoses how much water is retained in the foodplains and thus return to soil
711	!! moisture or is taken out of the rivers for irrigation. \n
712	!!
713	!! RECENT CHANGE(S): None
714	!!
715	!! MAIN OUTPUT VARIABLE(S):
716	!! The result of the routing are 3 fluxes :
717	!! - riverflow : The water which flows out from the major rivers. The flux will be located
718	!! on the continental grid but this should be a coastal point.
719	!! - coastalflow : This is the water which flows in a disperse way into the ocean. Essentially these
720	!! are the outflows from all of the small rivers.
721	!! - returnflow : This is the water which flows into a land-point - typically rivers which end in
722	!! the desert. This water will go back into the hydrol module to allow re-evaporation.
723	!! - irrigation : This is water taken from the reservoir and is being put into the upper
724	!! layers of the soil.
725	!! The two first fluxes are in kg/dt and the last two fluxes are in kg/(m^2dt).\n
726	!!
727	!! REFERENCE(S) :
728	!! - Miller JR, Russell GL, Caliri G (1994)
729	!! Continental-scale river flow in climate models.
730	!! J. Clim., 7:914-928
731	!! - Hagemann S and Dumenil L. (1998)
732	!! A parametrization of the lateral waterflow for the global scale.
733	!! Clim. Dyn., 14:17-31
734	!! - Oki, T., T. Nishimura, and P. Dirmeyer (1999)
735	!! Assessment of annual runoff from land surface models using total runoff integrating pathways (TRIP)
736	!! J. Meteorol. Soc. Jpn., 77, 235-255
737	!! - Fekete BM, Charles V, Grabs W (2000)
738	!! Global, composite runoff fields based on observed river discharge and simulated water balances.
739	!! Technical report, UNH/GRDC, Global Runoff Data Centre, Koblenz
740	!! - Vorosmarty, C., B. Fekete, B. Meybeck, and R. Lammers (2000)
741	!! Global system of rivers: Its role in organizing continental land mass and defining land-to-ocean linkages
742	!! Global Biogeochem. Cycles, 14, 599-621
743	!! - Vivant, A-C. (?? 2002)
744	!! DÃ©veloppement du schÃ©ma de routage et des plaines d'inondation, MSc Thesis, Paris VI University
745	!! - J. Polcher (2003)
746	!! Les processus de surface a l'echelle globale et leurs interactions avec l'atmosphere
747	!! Habilitation a diriger les recherches, Paris VI University, 67pp.
748	!!
749	!! FLOWCHART :
750	!! \latexonly
751	!! \includegraphics[scale=0.75]{routing_main_flowchart.png}
752	!! \endlatexonly
753	!! \n
754	!_ ================================================================================================================================
755
756	SUBROUTINE routing_highres_main(kjit, nbpt, index, &
757	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
758	& drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
759	& tempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
760
761	IMPLICIT NONE
762
763	!! 0.1 Input variables
764	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
765	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
766	INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless)
767	INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless)
768	INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless)
769	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless)
770	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
771	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb) !! Vector of neighbours for each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
772	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m)
773	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
774	REAL(r_std), INTENT(in) :: totfrac_nobio(nbpt) !! Total fraction of no-vegetation (continental ice, lakes ...) (unitless;0-1)
775	REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) !! Maximal fraction of vegetation (unitless;0-1)
776	REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt)
777	REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff (kg/m^2/dt)
778	REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt)
779	REAL(r_std), INTENT(in) :: transpot(nbpt,nvm) !! Potential transpiration of the vegetation (kg/m^2/dt)
780	REAL(r_std), INTENT(in) :: precip_rain(nbpt) !! Rainfall (kg/m^2/dt)
781	REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
782	REAL(r_std), INTENT(in) :: humrel(nbpt,nvm) !! Soil moisture stress, root extraction potential (unitless)
783	REAL(r_std), INTENT(in) :: tempdiag(nbpt,ngrnd) !! Diagnostic soil temperature profile
784	REAL(r_std), INTENT(in) :: reinf_slope(nbpt) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
785
786	!! 0.2 Output variables
787	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box.
788	!! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt)
789	REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
790	REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
791	REAL(r_std), INTENT(out) :: riverflow(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
792	REAL(r_std), INTENT(out) :: coastalflow(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
793	REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
794	REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
795
796	!! 0.3 Local variables
797	REAL(r_std), DIMENSION(nbpt) :: return_lakes !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
798
799	INTEGER(i_std) :: ig, jv !! Indices (unitless)
800	REAL(r_std), DIMENSION(nbpt) :: tot_vegfrac_nowoody !! Total fraction occupied by grass (0-1,unitless)
801
802	REAL(r_std), DIMENSION(nbpt) :: fast_diag_old !! Reservoir in the beginning of the time step
803	REAL(r_std), DIMENSION(nbpt) :: slow_diag_old !! Reservoir in the beginning of the time step
804	REAL(r_std), DIMENSION(nbpt) :: stream_diag_old !! Reservoir in the beginning of the time step
805	REAL(r_std), DIMENSION(nbpt) :: lake_diag_old !! Reservoir in the beginning of the time step
806	REAL(r_std), DIMENSION(nbpt) :: pond_diag_old !! Reservoir in the beginning of the time step
807	REAL(r_std), DIMENSION(nbpt) :: flood_diag_old !! Reservoir in the beginning of the time step
808
809	!! For water budget check in the three routing reservoirs (positive if input > output)
810	!! Net fluxes averaged over each grid cell in kg/m^2/dt
811	REAL(r_std), DIMENSION(nbpt) :: netflow_stream_diag !! Input - Output flow to stream reservoir
812	REAL(r_std), DIMENSION(nbpt) :: netflow_fast_diag !! Input - Output flow to fast reservoir
813	REAL(r_std), DIMENSION(nbpt) :: netflow_slow_diag !! Input - Output flow to slow reservoir
814	!
815	REAL(r_std), DIMENSION(nbpt,nbasmax) :: stemp_total_tend, stemp_advec_tend, stemp_relax_tend
816	!
817	!_ ================================================================================================================================
818
819	! Save reservoirs in beginning of time step to calculate the water budget
820	fast_diag_old = fast_diag
821	slow_diag_old = slow_diag
822	stream_diag_old = stream_diag
823	lake_diag_old = lake_diag
824	pond_diag_old = pond_diag
825	flood_diag_old = flood_diag
826
827	!
828	!! Computes the variables averaged between routing time steps and which will be used in subsequent calculations
829	!
830	floodout_mean(:) = floodout_mean(:) + floodout(:)
831	runoff_mean(:) = runoff_mean(:) + runoff(:)
832	drainage_mean(:) = drainage_mean(:) + drainage(:)
833	floodtemp(:) = tempdiag(:,floodtemp_lev)
834	precip_mean(:) = precip_mean(:) + precip_rain(:)
835	!
836	!! Computes the total fraction occupied by the grasses and the crops for each grid cell
837	tot_vegfrac_nowoody(:) = zero
838	DO jv = 1, nvm
839	IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
840	tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv)
841	END IF
842	END DO
843
844	DO ig = 1, nbpt
845	IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN
846	DO jv = 1,nvm
847	IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
848	transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig)
849	END IF
850	END DO
851	ELSE
852	IF (MAXVAL(veget_max(ig,2:nvm)) .GT. min_sechiba) THEN
853	DO jv = 2, nvm
854	transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,2:nvm))
855	ENDDO
856	ENDIF
857	ENDIF
858	ENDDO
859
860	!
861	! Averaged variables (i.e. *dt_sechiba/dt_routing). This accounts for the difference between the shorter
862	! timestep dt_sechiba of other parts of the model and the long dt_routing timestep (set to one day at present)
863	!
864	totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dt_sechiba/dt_routing
865	k_litt_mean(:) = k_litt_mean(:) + k_litt(:)*dt_sechiba/dt_routing
866	tempdiag_mean(:,:) = tempdiag_mean(:,:) + tempdiag(:,:)*dt_sechiba/dt_routing
867	!
868	! Only potentially vegetated surfaces are taken into account. At the start of
869	! the growing seasons we will give more weight to these areas.
870	!
871	DO jv=2,nvm
872	DO ig=1,nbpt
873	humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)veget_max(ig,jv)dt_sechiba/dt_routing
874	vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
875	ENDDO
876	ENDDO
877	!
878	time_counter = time_counter + dt_sechiba
879	!
880	! If the time has come we do the routing.
881	!
882	IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN
883	!
884	!! Computes the transport of water in the various reservoirs
885	!
886	CALL routing_hr_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, &
887	& vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, &
888	& tempdiag_mean, reinf_slope, lakeinflow_mean, returnflow_mean, reinfiltration_mean, &
889	& irrigation_mean, riverflow_mean, coastalflow_mean, hydrographs, slowflow_diag, flood_frac, &
890	& flood_res, netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, &
891	& stemp_total_tend, stemp_advec_tend, stemp_relax_tend)
892	!
893	!! Responsible for storing the water in lakes
894	!
895	CALL routing_hr_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, return_lakes)
896	!
897	returnflow_mean(:) = returnflow_mean(:) + return_lakes(:)
898
899	time_counter = zero
900	!
901	floodout_mean(:) = zero
902	runoff_mean(:) = zero
903	drainage_mean(:) = zero
904	transpot_mean(:) = zero
905	precip_mean(:) = zero
906	!
907	humrel_mean(:) = zero
908	totnobio_mean(:) = zero
909	k_litt_mean(:) = zero
910	tempdiag_mean(:,:) = zero
911	vegtot_mean(:) = zero
912
913	! Change the units of the routing fluxes from kg/dt_routing into kg/dt_sechiba
914	hydrographs(:) = hydrographs(:)/dt_routing*dt_sechiba
915	HTUhgmon(:,:) = HTUhgmon(:,:)/dt_routing*dt_sechiba
916	slowflow_diag(:) = slowflow_diag(:)/dt_routing*dt_sechiba
917
918	! Change the units of the routing fluxes from kg/m^2/dt_routing into kg/m^2/dt_sechiba
919	returnflow_mean(:) = returnflow_mean(:)/dt_routing*dt_sechiba
920	reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba
921	irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba
922	irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba
923
924	! Change units as above but at the same time transform the kg/dt_routing to m^3/dt_sechiba
925	riverflow_mean(:) = riverflow_mean(:)/dt_routing*dt_sechiba/mille
926	coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dt_sechiba/mille
927
928	! Water budget residu of the three routing reservoirs (in kg/m^2/s)
929	! Note that these diagnostics are done using local variables only calculated
930	! during the time steps when the routing is calculated
931	CALL xios_orchidee_send_field("wbr_stream",(stream_diag - stream_diag_old - netflow_stream_diag)/dt_routing)
932	CALL xios_orchidee_send_field("wbr_fast", (fast_diag - fast_diag_old - netflow_fast_diag)/dt_routing)
933	CALL xios_orchidee_send_field("wbr_slow", (slow_diag - slow_diag_old - netflow_slow_diag)/dt_routing)
934	CALL xios_orchidee_send_field("wbr_lake", (lake_diag - lake_diag_old - &
935	lakeinflow_mean + return_lakes)/dt_routing)
936	CALL xios_orchidee_send_field("StreamT_TotTend", stemp_total_tend)
937	CALL xios_orchidee_send_field("StreamT_AdvTend", stemp_advec_tend)
938	CALL xios_orchidee_send_field("StreamT_RelTend", stemp_relax_tend)
939	ENDIF
940
941	!
942	! Return the fraction of routed water for this time step.
943	!
944	returnflow(:) = returnflow_mean(:)
945	reinfiltration(:) = reinfiltration_mean(:)
946	irrigation(:) = irrigation_mean(:)
947	riverflow(:) = riverflow_mean(:)
948	coastalflow(:) = coastalflow_mean(:)
949
950	!
951	! Write diagnostics
952	!
953	!
954	CALL xios_orchidee_send_field("mask_coast",mask_coast)
955
956	IF ( do_irrigation ) THEN
957	CALL xios_orchidee_send_field("irrigmap",irrigated)
958	ENDIF
959
960	IF ( do_floodplains ) THEN
961	!! May be improved by performing the operation with XIOS
962	floodmap(:) = 0.0
963	DO ig=1,nbpt
964	floodmap(ig) = SUM(floodplains(ig,:)) / (area(ig)*contfrac(ig))
965	END DO
966	CALL xios_orchidee_send_field("floodmap",floodmap)
967	ENDIF
968
969	IF ( doswamps ) THEN
970	CALL xios_orchidee_send_field("swampmap",swamp)
971	ENDIF
972
973	!
974	! Water storage in reservoirs [kg/m^2]
975	CALL xios_orchidee_send_field("fastr",fast_diag)
976	CALL xios_orchidee_send_field("slowr",slow_diag)
977	CALL xios_orchidee_send_field("streamr",stream_diag)
978	CALL xios_orchidee_send_field("laker",lake_diag)
979	CALL xios_orchidee_send_field("pondr",pond_diag)
980	CALL xios_orchidee_send_field("floodr",flood_diag)
981	CALL xios_orchidee_send_field("floodh",flood_height)
982
983	! FLOODPLAINS
984	CALL xios_orchidee_send_field("flood_frac",flood_frac)
985
986	! Difference between the end and the beginning of the routing time step [kg/m^2]
987	CALL xios_orchidee_send_field("delfastr", fast_diag - fast_diag_old)
988	CALL xios_orchidee_send_field("delslowr", slow_diag - slow_diag_old)
989	CALL xios_orchidee_send_field("delstreamr", stream_diag - stream_diag_old)
990	CALL xios_orchidee_send_field("dellaker", lake_diag - lake_diag_old)
991	CALL xios_orchidee_send_field("delpondr", pond_diag - pond_diag_old)
992	CALL xios_orchidee_send_field("delfloodr", flood_diag - flood_diag_old)
993
994	! Water fluxes converted from kg/m^2/dt_sechiba into kg/m^2/s
995	CALL xios_orchidee_send_field("irrigation",irrigation/dt_sechiba)
996	CALL xios_orchidee_send_field("netirrig",irrig_netereq/dt_sechiba)
997	CALL xios_orchidee_send_field("riversret",returnflow/dt_sechiba)
998	CALL xios_orchidee_send_field("reinfiltration",reinfiltration/dt_sechiba)
999
1000	! Transform from kg/dt_sechiba into m^3/s
1001	CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba)
1002	CALL xios_orchidee_send_field("htuhgmon",HTUhgmon/mille/dt_sechiba)
1003	CALL xios_orchidee_send_field("htutempmon",HTUtempmon)
1004	CALL xios_orchidee_send_field("hydrotemp", hydrotemp)
1005	CALL xios_orchidee_send_field("streamlimit", streamlimit)
1006
1007	CALL xios_orchidee_send_field("slowflow",slowflow_diag/mille/dt_sechiba) ! previous id name: Qb
1008	CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
1009	CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
1010
1011	IF ( .NOT. xios_orchidee_ok) THEN
1012	IF ( .NOT. almaoutput ) THEN
1013	!
1014	CALL histwrite_p(hist_id, 'riversret', kjit, returnflow, nbpt, index)
1015	IF (do_floodplains .OR. doponds) THEN
1016	CALL histwrite_p(hist_id, 'reinfiltration', kjit, reinfiltration, nbpt, index)
1017	ENDIF
1018	CALL histwrite_p(hist_id, 'hydrographs', kjit, hydrographs/mille, nbpt, index)
1019	!
1020	CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index)
1021	CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index)
1022	CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index)
1023	IF ( do_floodplains ) THEN
1024	CALL histwrite_p(hist_id, 'floodr', kjit, flood_diag, nbpt, index)
1025	CALL histwrite_p(hist_id, 'floodh', kjit, flood_height, nbpt, index)
1026	ENDIF
1027	CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, nbpt, index)
1028	CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index)
1029	!
1030	IF ( do_irrigation ) THEN
1031	CALL histwrite_p(hist_id, 'irrigation', kjit, irrigation, nbpt, index)
1032	CALL histwrite_p(hist_id, 'returnflow', kjit, returnflow, nbpt, index)
1033	CALL histwrite_p(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
1034	ENDIF
1035	!
1036	ELSE
1037	CALL histwrite_p(hist_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index)
1038	CALL histwrite_p(hist_id, 'Dis', kjit, hydrographs/mille, nbpt, index)
1039	!
1040	CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index)
1041	CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index)
1042	CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index)
1043	CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index)
1044	CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, nbpt, index)
1045	!
1046	IF ( do_irrigation ) THEN
1047	CALL histwrite_p(hist_id, 'Qirrig', kjit, irrigation, nbpt, index)
1048	CALL histwrite_p(hist_id, 'Qirrig_req', kjit, irrig_netereq, nbpt, index)
1049	ENDIF
1050	!
1051	ENDIF
1052	IF ( hist2_id > 0 ) THEN
1053	IF ( .NOT. almaoutput) THEN
1054	!
1055	CALL histwrite_p(hist2_id, 'riversret', kjit, returnflow, nbpt, index)
1056	IF (do_floodplains .OR. doponds) THEN
1057	CALL histwrite_p(hist2_id, 'reinfiltration', kjit, reinfiltration, nbpt, index)
1058	ENDIF
1059	CALL histwrite_p(hist2_id, 'hydrographs', kjit, hydrographs/mille, nbpt, index)
1060	!
1061	CALL histwrite_p(hist2_id, 'fastr', kjit, fast_diag, nbpt, index)
1062	CALL histwrite_p(hist2_id, 'slowr', kjit, slow_diag, nbpt, index)
1063	IF ( do_floodplains ) THEN
1064	CALL histwrite_p(hist2_id, 'floodr', kjit, flood_diag, nbpt, index)
1065	CALL histwrite_p(hist2_id, 'floodh', kjit, flood_height, nbpt, index)
1066	ENDIF
1067	CALL histwrite_p(hist2_id, 'pondr', kjit, pond_diag, nbpt, index)
1068	CALL histwrite_p(hist2_id, 'streamr', kjit, stream_diag, nbpt, index)
1069	CALL histwrite_p(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index)
1070	!
1071	IF ( do_irrigation ) THEN
1072	CALL histwrite_p(hist2_id, 'irrigation', kjit, irrigation, nbpt, index)
1073	CALL histwrite_p(hist2_id, 'returnflow', kjit, returnflow, nbpt, index)
1074	CALL histwrite_p(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
1075	ENDIF
1076	!
1077	ELSE
1078	!
1079	CALL histwrite_p(hist2_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index)
1080	CALL histwrite_p(hist2_id, 'Dis', kjit, hydrographs/mille, nbpt, index)
1081	!
1082	ENDIF
1083	ENDIF
1084	ENDIF
1085	!
1086	!
1087	END SUBROUTINE routing_highres_main
1088
1089	!! =============================================================================================================================
1090	!! SUBROUTINE: routing_highres_finalize
1091	!!
1092	!>\BRIEF Write to restart file
1093	!!
1094	!! DESCRIPTION: Write module variables to restart file
1095	!!
1096	!! RECENT CHANGE(S)
1097	!!
1098	!! REFERENCE(S)
1099	!!
1100	!! FLOWCHART
1101	!! \n
1102	!_ ==============================================================================================================================
1103
1104	SUBROUTINE routing_highres_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
1105
1106	IMPLICIT NONE
1107
1108	!! 0.1 Input variables
1109	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1110	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
1111	INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless)
1112	REAL(r_std), INTENT(in) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
1113	REAL(r_std), INTENT(in) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
1114
1115	!! 0.2 Local variables
1116
1117	!_ ================================================================================================================================
1118
1119	!
1120	! Write restart variables
1121	!
1122	CALL restput_p (rest_id, 'routingcounter', kjit, time_counter)
1123
1124	CALL restput_p (rest_id, 'streamtcst', kjit, stream_tcst)
1125	CALL restput_p (rest_id, 'slowtcst', kjit, slow_tcst)
1126	CALL restput_p (rest_id, 'fasttcst', kjit, fast_tcst)
1127	CALL restput_p (rest_id, 'floodtcst', kjit, flood_tcst)
1128	CALL restput_p (rest_id, 'swampcst', kjit, swamp_cst)
1129
1130	CALL restput_p (rest_id, 'lim_floodcri', kjit, lim_floodcri)
1131
1132	CALL restput_p (rest_id, 'nbasmax', kjit, nbasmax)
1133	CALL restput_p (rest_id, 'nbasmon', kjit, nbasmon)
1134	CALL restput_p (rest_id, 'inflows', kjit, inflows)
1135
1136	CALL restput_p (rest_id, 'routingarea', nbp_glo, nbasmax, 1, kjit, routing_area, 'scatter', nbp_glo, index_g)
1137	CALL restput_p (rest_id, 'routetogrid', nbp_glo, nbasmax, 1, kjit, REAL(route_togrid,r_std), 'scatter', &
1138	nbp_glo, index_g)
1139	CALL restput_p (rest_id, 'routetobasin', nbp_glo, nbasmax, 1, kjit, REAL(route_tobasin,r_std), 'scatter', &
1140	nbp_glo, index_g)
1141	CALL restput_p (rest_id, 'routenbintobas', nbp_glo, nbasmax, 1, kjit, REAL(route_nbintobas,r_std), 'scatter', &
1142	nbp_glo, index_g)
1143	CALL restput_p (rest_id, 'basinid', nbp_glo, nbasmax, 1, kjit, REAL(global_basinid,r_std), 'scatter', &
1144	nbp_glo, index_g)
1145	CALL restput_p (rest_id, 'topoindex', nbp_glo, nbasmax, 1, kjit, topo_resid, 'scatter', nbp_glo, index_g)
1146	CALL restput_p (rest_id, 'topoindex_stream', nbp_glo, nbasmax, 1, kjit, stream_resid, 'scatter', nbp_glo, index_g)
1147	CALL restput_p (rest_id, 'fastres', nbp_glo, nbasmax, 1, kjit, fast_reservoir, 'scatter', nbp_glo, index_g)
1148	CALL restput_p (rest_id, 'slowres', nbp_glo, nbasmax, 1, kjit, slow_reservoir, 'scatter', nbp_glo, index_g)
1149	CALL restput_p (rest_id, 'streamres', nbp_glo, nbasmax, 1, kjit, stream_reservoir, 'scatter',nbp_glo,index_g)
1150	CALL restput_p (rest_id, 'floodres', nbp_glo, nbasmax, 1, kjit, flood_reservoir, 'scatter', nbp_glo, index_g)
1151	CALL restput_p (rest_id, 'floodh', nbp_glo, nbasmax, 1, kjit, flood_height, 'scatter', nbp_glo, index_g)
1152	CALL restput_p (rest_id, 'flood_frac_bas', nbp_glo, nbasmax, 1, kjit, flood_frac_bas, 'scatter', nbp_glo, index_g)
1153	CALL restput_p (rest_id, 'pond_frac', nbp_glo, 1, 1, kjit, pond_frac, 'scatter', nbp_glo, index_g)
1154	CALL restput_p (rest_id, 'flood_frac', nbp_glo, 1, 1, kjit, flood_frac, 'scatter', nbp_glo, index_g)
1155	CALL restput_p (rest_id, 'flood_res', nbp_glo, 1, 1, kjit, flood_res, 'scatter', nbp_glo, index_g)
1156
1157	CALL restput_p (rest_id, 'fasttemp', nbp_glo, nbasmax, 1, kjit, fast_temp, 'scatter', nbp_glo, index_g)
1158	CALL restput_p (rest_id, 'slowtemp', nbp_glo, nbasmax, 1, kjit, slow_temp, 'scatter', nbp_glo, index_g)
1159	CALL restput_p (rest_id, 'streamtemp', nbp_glo, nbasmax, 1, kjit, stream_temp, 'scatter',nbp_glo,index_g)
1160
1161
1162	CALL restput_p (rest_id, 'lakeres', nbp_glo, 1, 1, kjit, lake_reservoir, 'scatter', nbp_glo, index_g)
1163	CALL restput_p (rest_id, 'pondres', nbp_glo, 1, 1, kjit, pond_reservoir, 'scatter', nbp_glo, index_g)
1164
1165	CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g)
1166	CALL restput_p (rest_id, 'returnflow', nbp_glo, 1, 1, kjit, returnflow_mean, 'scatter', nbp_glo, index_g)
1167	CALL restput_p (rest_id, 'reinfiltration', nbp_glo, 1, 1, kjit, reinfiltration_mean, 'scatter', nbp_glo, index_g)
1168	CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g)
1169	CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g)
1170	CALL restput_p (rest_id, 'hydrographs', nbp_glo, 1, 1, kjit, hydrographs, 'scatter', nbp_glo, index_g)
1171	CALL restput_p (rest_id, 'htuhgmon', nbp_glo, nbasmon, 1, kjit, HTUhgmon, 'scatter', nbp_glo, index_g)
1172	CALL restput_p (rest_id, 'slowflow_diag', nbp_glo, 1, 1, kjit, slowflow_diag, 'scatter', nbp_glo, index_g)
1173	CALL restput_p (rest_id, 'hydrotemp', nbp_glo, 1, 1, kjit, hydrotemp, 'scatter', nbp_glo, index_g)
1174	CALL restput_p (rest_id, 'htutempmon', nbp_glo, nbasmon, 1, kjit, HTUtempmon, 'scatter', nbp_glo, index_g)
1175	!
1176	! Keep track of the accumulated variables
1177	!
1178	CALL restput_p (rest_id, 'floodout_route', nbp_glo, 1, 1, kjit, floodout_mean, 'scatter', nbp_glo, index_g)
1179	CALL restput_p (rest_id, 'runoff_route', nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g)
1180	CALL restput_p (rest_id, 'drainage_route', nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g)
1181	CALL restput_p (rest_id, 'transpot_route', nbp_glo, 1, 1, kjit, transpot_mean, 'scatter', nbp_glo, index_g)
1182	CALL restput_p (rest_id, 'precip_route', nbp_glo, 1, 1, kjit, precip_mean, 'scatter', nbp_glo, index_g)
1183	CALL restput_p (rest_id, 'humrel_route', nbp_glo, 1, 1, kjit, humrel_mean, 'scatter', nbp_glo, index_g)
1184	CALL restput_p (rest_id, 'totnobio_route', nbp_glo, 1, 1, kjit, totnobio_mean, 'scatter', nbp_glo, index_g)
1185	CALL restput_p (rest_id, 'k_litt_route', nbp_glo, 1, 1, kjit, k_litt_mean, 'scatter', nbp_glo, index_g)
1186	CALL restput_p (rest_id, 'vegtot_route', nbp_glo, 1, 1, kjit, vegtot_mean, 'scatter', nbp_glo, index_g)
1187	CALL restput_p (rest_id, 'tempdiag_route', nbp_glo, ngrnd, 1, kjit, tempdiag_mean, 'scatter', nbp_glo, index_g)
1188
1189	CALL restput_p (rest_id, 'gridrephtu', nbp_glo, 1, 1, kjit, REAL(hydrodiag,r_std), 'scatter', nbp_glo, index_g)
1190	CALL restput_p (rest_id, 'htudiag', nbp_glo, nbasmon, 1, kjit, REAL(HTUdiag_loc,r_std), 'scatter', nbp_glo, index_g)
1191
1192	IF ( do_irrigation ) THEN
1193	CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter', nbp_glo, index_g)
1194	CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter', nbp_glo, index_g)
1195	ENDIF
1196
1197	IF ( do_floodplains ) THEN
1198	CALL restput_p (rest_id, 'floodplains', nbp_glo, nbasmax, 1, kjit, floodplains, 'scatter', nbp_glo, index_g)
1199	CALL restput_p (rest_id, 'floodcri', nbp_glo, nbasmax, 1, kjit, floodcri, 'scatter', nbp_glo, index_g)
1200	CALL restput_p (rest_id, 'floodp_beta', nbp_glo, nbasmax, 1, kjit, fp_beta, 'scatter', nbp_glo, index_g)
1201	ENDIF
1202	IF ( dofloodoverflow ) THEN
1203	CALL restput_p (rest_id, 'orog_min', nbp_glo, nbasmax, 1,kjit,orog_min, 'scatter', nbp_glo, index_g)
1204	END IF
1205	IF ( doswamps ) THEN
1206	CALL restput_p (rest_id, 'swamp', nbp_glo, 1, 1, kjit, swamp, 'scatter', nbp_glo, index_g)
1207	ENDIF
1208
1209	END SUBROUTINE routing_highres_finalize
1210
1211	!! ================================================================================================================================
1212	!! SUBROUTINE : routing_hr_init
1213	!!
1214	!>\BRIEF This subroutine allocates the memory and get the fixed fields from the restart file.
1215	!!
1216	!! DESCRIPTION (definitions, functional, design, flags) : None
1217	!!
1218	!! RECENT CHANGE(S): None
1219	!!
1220	!! MAIN OUTPUT VARIABLE(S):
1221	!!
1222	!! REFERENCES : None
1223	!!
1224	!! FLOWCHART :None
1225	!! \n
1226	!_ ================================================================================================================================
1227
1228	SUBROUTINE routing_hr_init(kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
1229	& riverflow, coastalflow, flood_frac, flood_res, tempdiag, rest_id)
1230	!
1231	IMPLICIT NONE
1232	!
1233	! interface description
1234	!
1235	!! INPUT VARIABLES
1236	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1237	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
1238	INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index !! Indices of the points on the map (unitless)
1239	REAL(r_std), DIMENSION(nbpt,ngrnd),INTENT(in) :: tempdiag !! Temperature profile in soil
1240	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless)
1241	!
1242	!! OUTPUT VARIABLES
1243	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: returnflow !! The water flow from lakes and swamps which returns into the grid box.
1244	!! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt)
1245	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: reinfiltration !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
1246	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: irrigation !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil.(kg/m^2/dt)
1247	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: riverflow !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
1248	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: coastalflow !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
1249	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_frac !! Flooded fraction of the grid box (unitless;0-1)
1250	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_res !! Diagnostic of water amount in the floodplains reservoir (kg)
1251	!
1252	!! LOCAL VARIABLES
1253	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
1254	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g !! A temporary real array for the integers
1255	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: tmp_real !
1256	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g2
1257	REAL(r_std) :: ratio !! Diagnostic ratio to check that dt_routing is a multiple of dt_sechiba (unitless)
1258	REAL(r_std) :: totarea !! Total area of basin (m^2)
1259	INTEGER(i_std) :: ier, ig, im, ib, ipn(1), nbhtumon !! Indices (unitless)
1260	REAL(r_std) :: nbasmon_tmp, nbasmax_tmp, inflows_tmp
1261
1262	!_ ================================================================================================================================
1263	!
1264	!
1265	! These variables will require the configuration infrastructure
1266	!
1267	!Config Key = DT_ROUTING
1268	!Config If = RIVER_ROUTING
1269	!Config Desc = Time step of the routing scheme
1270	!Config Def = one_day
1271	!Config Help = This values gives the time step in seconds of the routing scheme.
1272	!Config It should be multiple of the main time step of ORCHIDEE. One day
1273	!Config is a good value.
1274	!Config Units = [seconds]
1275	!
1276	dt_routing = dt_sechiba
1277	CALL getin_p('DT_ROUTING', dt_routing)
1278	!
1279	!
1280	!
1281	!Config Key = DO_FLOODINFILT
1282	!Config Desc = Should floodplains reinfiltrate into the soil
1283	!Config If = RIVER_ROUTING
1284	!Config Def = n
1285	!Config Help = This parameters allows the user to ask the model
1286	!Config to take into account the flood plains reinfiltration
1287	!Config into the soil moisture. It then can go
1288	!Config back to the slow and fast reservoirs
1289	!Config Units = [FLAG]
1290	!
1291	dofloodinfilt = .FALSE.
1292	IF ( do_floodplains ) CALL getin_p('DO_FLOODINFILT', dofloodinfilt)
1293	!
1294	!Config Key = CONDUCT_FACTOR_FP
1295	!Config Desc = Adjustment factor for floodplains reinfiltration
1296	!Config If = RIVER_ROUTING
1297	!Config Def = n
1298	!Config Help = Factor used to reduce the infiltration from the
1299	!Config floodplains. For a value of 1, the infiltration is
1300	!Config unchanged, for a value of 0 there is no infiltration.
1301	!Config Units = -
1302	!
1303	conduct_factor = 1.0
1304	IF ( do_floodplains ) CALL getin_p('CONDUCT_FACTOR_FP', conduct_factor)
1305	!
1306	!
1307	!Config Key = DO_FLOODOVERFLOW
1308	!Config Desc = Should floodplains overflow to upstream HTUs floodplains
1309	!Config If = RIVER_ROUTING
1310	!Config Def = n
1311	!Config Help = This parameters allows the user to ask the model
1312	!Config to take into account the overflow of the
1313	!Config floodplains. The water can flow to the upstream
1314	!Config floodplains reservoir if the current flood height
1315	!Config is higher than the upstream one.
1316	!Config Units = [FLAG]
1317	!
1318	dofloodoverflow = .FALSE.
1319	IF ( do_floodplains ) CALL getin_p('DO_FLOODOVERFLOW', dofloodoverflow)
1320	!
1321	!Config Key = OVERFLOW_REPETITION
1322	!Config Desc = Repetition of overflow at each routing time step
1323	!Config If = RIVER_ROUTING
1324	!Config Def = n
1325	!Config Help = This parameters allows the user to ask the model
1326	!Config repeat the overflow a certain amount of time
1327	!Config in order to have more stability with lower
1328	!Config overflow time step.
1329	!Config Units = [FLAG]
1330	!
1331	overflow_repetition = 1
1332	IF ( do_floodplains ) CALL getin_p('OVERFLOW_REPETITION', overflow_repetition)
1333	!
1334	!Config Key = R_FLOODMAX
1335	!Config Desc = Maximal values for R factor
1336	!Config If = DO_FLOODPLAINS
1337	!Config Def = 0.5
1338	!Config Help = R is the factor of reduction of the stream discharge
1339	!Config if there is floodplains. This is the maximal value
1340	!Config when the HTU is fully filled.
1341	!Config R = 1 -> discharge = 0
1342	!Config R = 0 -> Maximal discharge
1343	!
1344	rfloodmax = 0.5
1345	IF ( do_floodplains ) CALL getin_p('R_FLOODMAX', rfloodmax)
1346	!
1347	!Config Key = OVERFLOW_TCST
1348	!Config Desc = Time Constant for overflow in day
1349	!Config If = DO_FLOODPLAINS
1350	!Config Def = 1
1351	!Config Help = OVERFLOW_TCST is the time constant
1352	!Config For the floodplains overflow
1353	!
1354	overflow_tcst = 1
1355	IF ( do_floodplains ) CALL getin_p('OVERFLOW_TCST', overflow_tcst)
1356	!
1357	!Config Key = DO_SWAMPS
1358	!Config Desc = Should we include swamp parameterization
1359	!Config If = RIVER_ROUTING
1360	!Config Def = n
1361	!Config Help = This parameters allows the user to ask the model
1362	!Config to take into account the swamps and return
1363	!Config the water into the bottom of the soil. It then can go
1364	!Config back to the atmopshere. This tried to simulate
1365	!Config internal deltas of rivers.
1366	!Config Units = [FLAG]
1367	!
1368	doswamps = .FALSE.
1369	CALL getin_p('DO_SWAMPS', doswamps)
1370	!
1371	!Config Key = DO_PONDS
1372	!Config Desc = Should we include ponds
1373	!Config If = RIVER_ROUTING
1374	!Config Def = n
1375	!Config Help = This parameters allows the user to ask the model
1376	!Config to take into account the ponds and return
1377	!Config the water into the soil moisture. It then can go
1378	!Config back to the atmopshere. This tried to simulate
1379	!Config little ponds especially in West Africa.
1380	!Config Units = [FLAG]
1381	!
1382	doponds = .FALSE.
1383	CALL getin_p('DO_PONDS', doponds)
1384	!
1385	!Config Key = FLOOD_BETA
1386	!Config Desc = Parameter to fix the shape of the floodplain
1387	!Config If = RIVER_ROUTING
1388	!Config Def = 2.0
1389	!Config Help = Parameter to fix the shape of the floodplain
1390	!Config (>1 for convex edges, <1 for concave edges)
1391	!Config Units = [-]
1392
1393	! ANTHONY OLD FLOODPLAINS
1394	!CALL getin_p("FLOOD_BETA", beta)
1395	!
1396	!Config Key = POND_BETAP
1397	!Config Desc = Ratio of the basin surface intercepted by ponds and the maximum surface of ponds
1398	!Config If = RIVER_ROUTING
1399	!Config Def = 0.5
1400	!Config Help =
1401	!Config Units = [-]
1402	CALL getin_p("POND_BETAP", betap)
1403	!
1404	!Config Key = FLOOD_CRI
1405	!Config Desc = Potential height for which all the basin is flooded
1406	!Config If = DO_FLOODPLAINS or DO_PONDS
1407	!Config Def = 2000.
1408	!Config Help =
1409	!Config Units = [mm]
1410
1411	! ANTHONY OLD FLOODPLAINS
1412	!CALL getin_p("FLOOD_CRI", floodcri)
1413	!
1414	!Config Key = POND_CRI
1415	!Config Desc = Potential height for which all the basin is a pond
1416	!Config If = DO_FLOODPLAINS or DO_PONDS
1417	!Config Def = 2000.
1418	!Config Help =
1419	!Config Units = [mm]
1420	CALL getin_p("POND_CRI", pondcri)
1421
1422	!Config Key = MAX_LAKE_RESERVOIR
1423	!Config Desc = Maximum limit of water in lake_reservoir
1424	!Config If = RIVER_ROUTING
1425	!Config Def = 7000
1426	!Config Help =
1427	!Config Units = [kg/m2(routing area)]
1428	max_lake_reservoir = 7000
1429	CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir)
1430
1431	!
1432	!
1433	! In order to simplify the time cascade check that dt_routing
1434	! is a multiple of dt_sechiba
1435	!
1436	ratio = dt_routing/dt_sechiba
1437	IF ( ABS(NINT(ratio) - ratio) .GT. 10*EPSILON(ratio)) THEN
1438	WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
1439	WRITE(numout,*) "The chosen time step for the routing is not a multiple of the"
1440	WRITE(numout,*) "main time step of the model. We will change dt_routing so that"
1441	WRITE(numout,*) "this condition os fulfilled"
1442	dt_routing = NINT(ratio) * dt_sechiba
1443	WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
1444	ENDIF
1445	!
1446	IF ( dt_routing .LT. dt_sechiba) THEN
1447	WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
1448	WRITE(numout,*) 'The routing timestep can not be smaller than the one'
1449	WRITE(numout,*) 'of the model. We reset its value to the model''s timestep.'
1450	WRITE(numout,*) 'The old DT_ROUTING is : ', dt_routing
1451	dt_routing = dt_sechiba
1452	WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
1453	ENDIF
1454	!
1455	! If the routing_graph file is available we will extract the information in the dimensions
1456	! and parameters.
1457	!
1458	!Config Key = ROUTING_FILE
1459	!Config Desc = Name of file which contains the routing information graph on the model grid
1460	!Config If = RIVER_ROUTING
1461	!Config Def = routing.nc
1462	!Config Help = The file provided here should allows to route the water from one HTU
1463	!Config to another. The RoutingPP code needs to be used in order to generate
1464	!Config the routing graph for the model grid.
1465	!Config More details on : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp
1466	!Config Units = [FILE]
1467	!
1468	!graphfilename = 'routing_graph.nc'
1469	!CALL getin('ROUTING_FILE',graphfilename)
1470	!CALL routing_hr_graphinfo(graphfilename, nbasmax, inflows, nbasmon, undef_graphfile, stream_tcst, fast_tcst, slow_tcst, & & flood_tcst, swamp_cst, lim_floodcri)
1471	! At this stage we could have an option to force reading of graph
1472	!
1473	! Constants which can be in the restart file
1474	!
1475	var_name ="routingcounter"
1476	CALL ioconf_setatt_p('UNITS', 's')
1477	CALL ioconf_setatt_p('LONG_NAME','Time counter for the routing scheme')
1478	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, time_counter)
1479	CALL setvar_p (time_counter, val_exp, 'NO_KEYWORD', zero)
1480
1481	! Parameters which are in the restart file
1482	IF (stream_tcst .LE. 0 ) THEN
1483	var_name ="streamtcst"
1484	CALL ioconf_setatt_p('UNITS', 's/km')
1485	CALL ioconf_setatt_p('LONG_NAME','Time constant for the stream reservoir')
1486	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, stream_tcst)
1487	ENDIF
1488	IF (slow_tcst .LE. 0 ) THEN
1489	var_name ="slowtcst"
1490	CALL ioconf_setatt_p('UNITS', 's/km')
1491	CALL ioconf_setatt_p('LONG_NAME','Time constant for the slow reservoir')
1492	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, slow_tcst)
1493	ENDIF
1494	IF (fast_tcst .LE. 0 ) THEN
1495	var_name ="fasttcst"
1496	CALL ioconf_setatt_p('UNITS', 's/km')
1497	CALL ioconf_setatt_p('LONG_NAME','Time constant for the fast reservoir')
1498	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, fast_tcst)
1499	ENDIF
1500	IF (flood_tcst .LE. 0 ) THEN
1501	var_name ="floodtcst"
1502	CALL ioconf_setatt_p('UNITS', 's/km')
1503	CALL ioconf_setatt_p('LONG_NAME','Time constant for the flood reservoir')
1504	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, flood_tcst)
1505	ENDIF
1506	IF (swamp_cst .LE. 0 ) THEN
1507	var_name ="swampcst"
1508	CALL ioconf_setatt_p('UNITS', '-')
1509	CALL ioconf_setatt_p('LONG_NAME','Fraction of the river transport that flows to the swamps')
1510	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, swamp_cst)
1511	ENDIF
1512	IF (lim_floodcri .LE. 0 ) THEN
1513	var_name ="lim_floodcri"
1514	CALL ioconf_setatt_p('UNITS', 'm')
1515	CALL ioconf_setatt_p('LONG_NAME','Minimal difference of orography consecutive floodplains HTUs')
1516	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, lim_floodcri)
1517	ENDIF
1518	!
1519	! Number of HTUs
1520	!
1521	var_name ="nbasmax"
1522	CALL ioconf_setatt_p('UNITS', '-')
1523	CALL ioconf_setatt_p('LONG_NAME','Number of HTU per grid box')
1524	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, nbasmax_tmp)
1525	CALL routing_hr_restartconsistency(var_name, nbasmax, nbasmax_tmp)
1526	!
1527	! Number of inflows
1528	!
1529	var_name ="inflows"
1530	CALL ioconf_setatt_p('UNITS', '-')
1531	CALL ioconf_setatt_p('LONG_NAME','Maximum number of inflows per HTU')
1532	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, inflows_tmp)
1533	CALL routing_hr_restartconsistency(var_name, inflows, inflows_tmp)
1534	!
1535	! Dimension of HTU monitoring variable
1536	!
1537	var_name ="nbasmon"
1538	CALL ioconf_setatt_p('UNITS', '-')
1539	CALL ioconf_setatt_p('LONG_NAME','Number of HTU to be monitored')
1540	CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, nbasmon_tmp)
1541	CALL routing_hr_restartconsistency(var_name, nbasmon, nbasmon_tmp)
1542	!
1543	! Continuation of extraction from restart file.
1544	!
1545	ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier)
1546	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for routing_area_loc','','')
1547
1548	ALLOCATE (routing_area_glo(nbp_glo,nbasmax))
1549	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for routing_area_glo','','')
1550	IF ( .NOT. ReadGraph ) THEN
1551	var_name = 'routingarea'
1552	IF (is_root_prc) THEN
1553	CALL ioconf_setatt('UNITS', 'm^2')
1554	CALL ioconf_setatt('LONG_NAME','Area of basin')
1555	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g)
1556	ENDIF
1557	CALL scatter(routing_area_glo,routing_area_loc)
1558	routing_area=>routing_area_loc
1559	ENDIF
1560	CALL scatter(routing_area_glo,routing_area_loc)
1561	routing_area=>routing_area_loc
1562
1563
1564	IF ( do_floodplains ) THEN
1565	!!!!!!!!!!!!!!!!!!!!!!!!!!!
1566	!! ANTHONY - BETA
1567	ALLOCATE (fp_beta_loc(nbpt,nbasmax), stat=ier)
1568	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fp_beta_loc','','')
1569
1570	ALLOCATE (fp_beta_glo(nbp_glo,nbasmax))
1571	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fp_beta_glo','','')
1572
1573	IF ( .NOT. ReadGraph ) THEN
1574	IF (is_root_prc) THEN
1575	var_name = 'floodp_beta'
1576	CALL ioconf_setatt('UNITS', '-')
1577	CALL ioconf_setatt('LONG_NAME','Beta parameter for floodplains')
1578	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fp_beta_glo, "gather", nbp_glo, index_g)
1579	ENDIF
1580	CALL scatter(fp_beta_glo,fp_beta_loc)
1581	fp_beta=>fp_beta_loc
1582	END IF
1583	!!!!!!!!!!!!!!!!!!!!!!!!!!!
1584
1585	!!!!!!!!!!!!!!!!!!!!!!!!!!!
1586	!! ANTHONY - h0 - floodcri
1587	ALLOCATE (floodcri_loc(nbpt,nbasmax), stat=ier)
1588	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodcri_loc','','')
1589
1590	ALLOCATE (floodcri_glo(nbp_glo,nbasmax))
1591	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodcri_glo','','')
1592
1593	IF ( .NOT. ReadGraph ) THEN
1594	IF (is_root_prc) THEN
1595	var_name = 'floodcri'
1596	CALL ioconf_setatt('UNITS', 'mm')
1597	CALL ioconf_setatt('LONG_NAME','Height of complete flood')
1598	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., floodcri_glo, "gather", nbp_glo, index_g)
1599	END IF
1600	CALL scatter(floodcri_glo,floodcri_loc)
1601	floodcri=>floodcri_loc
1602	ENDIF
1603	END IF
1604
1605	ALLOCATE (tmp_real_g(nbp_glo,nbasmax), stat=ier)
1606	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tmp_real_g','','')
1607
1608	ALLOCATE (route_togrid_loc(nbpt,nbasmax), stat=ier)
1609	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_togrid_loc','','')
1610	ALLOCATE (route_togrid_glo(nbp_glo,nbasmax), stat=ier) ! used in global in routing_hr_flow
1611	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_togrid_glo','','')
1612
1613	IF ( .NOT. ReadGraph ) THEN
1614	IF (is_root_prc) THEN
1615	var_name = 'routetogrid'
1616	CALL ioconf_setatt('UNITS', '-')
1617	CALL ioconf_setatt('LONG_NAME','Grid into which the basin flows')
1618	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
1619	route_togrid_glo(:,:) = undef_int
1620	WHERE ( tmp_real_g .LT. val_exp )
1621	route_togrid_glo = NINT(tmp_real_g)
1622	ENDWHERE
1623	ENDIF
1624	CALL bcast(route_togrid_glo) ! used in global in routing_hr_flow
1625	CALL scatter(route_togrid_glo,route_togrid_loc)
1626	route_togrid=>route_togrid_loc
1627	ENDIF
1628	!
1629	ALLOCATE (route_tobasin_loc(nbpt,nbasmax), stat=ier)
1630	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_tobasin_loc','','')
1631
1632	ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax), stat=ier)
1633	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_tobasin_glo','','')
1634
1635	IF ( .NOT. ReadGraph ) THEN
1636	IF (is_root_prc) THEN
1637	var_name = 'routetobasin'
1638	CALL ioconf_setatt('UNITS', '-')
1639	CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes')
1640	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
1641	route_tobasin_glo = undef_int
1642	WHERE ( tmp_real_g .LT. val_exp )
1643	route_tobasin_glo = NINT(tmp_real_g)
1644	ENDWHERE
1645	num_largest = COUNT(route_tobasin_glo .EQ. nbasmax+3)
1646	ENDIF
1647	CALL scatter(route_tobasin_glo,route_tobasin_loc)
1648	CALL bcast(num_largest)
1649	route_tobasin=>route_tobasin_loc
1650	ENDIF
1651	!
1652	! nbintobasin
1653	!
1654	ALLOCATE (route_nbintobas_loc(nbpt,nbasmax), stat=ier)
1655	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_nbintobas_loc','','')
1656	ALLOCATE (route_nbintobas_glo(nbp_glo,nbasmax), stat=ier)
1657	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_nbintobas_glo','','')
1658
1659	IF ( .NOT. ReadGraph ) THEN
1660	IF (is_root_prc) THEN
1661	var_name = 'routenbintobas'
1662	CALL ioconf_setatt('UNITS', '-')
1663	CALL ioconf_setatt('LONG_NAME','Number of basin into current one')
1664	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
1665	route_nbintobas_glo = undef_int
1666	WHERE ( tmp_real_g .LT. val_exp )
1667	route_nbintobas_glo = NINT(tmp_real_g)
1668	ENDWHERE
1669	ENDIF
1670	CALL scatter(route_nbintobas_glo,route_nbintobas_loc)
1671	route_nbintobas=>route_nbintobas_loc
1672	ENDIF
1673	!
1674	ALLOCATE (global_basinid_loc(nbpt,nbasmax), stat=ier)
1675	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for global_basinid_loc','','')
1676	ALLOCATE (global_basinid_glo(nbp_glo,nbasmax), stat=ier)
1677	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for global_basinid_glo','','')
1678
1679	IF ( .NOT. ReadGraph ) THEN
1680	IF (is_root_prc) THEN
1681	var_name = 'basinid'
1682	CALL ioconf_setatt('UNITS', '-')
1683	CALL ioconf_setatt('LONG_NAME','ID of basin')
1684	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
1685	global_basinid_glo = undef_int
1686	WHERE ( tmp_real_g .LT. val_exp )
1687	global_basinid_glo = NINT(tmp_real_g)
1688	ENDWHERE
1689	ENDIF
1690	CALL scatter(global_basinid_glo,global_basinid_loc)
1691	global_basinid=>global_basinid_loc
1692	ENDIF
1693	!
1694	ALLOCATE (topo_resid_loc(nbpt,nbasmax), stat=ier)
1695	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for topo_resid_loc','','')
1696	ALLOCATE (topo_resid_glo(nbp_glo,nbasmax), stat=ier)
1697	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for topo_resid_glo','','')
1698
1699	IF ( .NOT. ReadGraph ) THEN
1700	IF (is_root_prc) THEN
1701	var_name = 'topoindex'
1702	CALL ioconf_setatt('UNITS', 'km')
1703	CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time')
1704	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g)
1705	ENDIF
1706	CALL scatter(topo_resid_glo,topo_resid_loc)
1707	topo_resid=>topo_resid_loc
1708	ENDIF
1709	!
1710	ALLOCATE (stream_resid_loc(nbpt,nbasmax), stat=ier)
1711	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_resid_loc','','')
1712	ALLOCATE (stream_resid_glo(nbp_glo,nbasmax), stat=ier)
1713	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_resid_glo','','')
1714
1715	IF ( .NOT. ReadGraph ) THEN
1716	IF (is_root_prc) THEN
1717	var_name = 'topoindex_stream'
1718	CALL ioconf_setatt('UNITS', 'km')
1719	CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time')
1720	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_resid_glo, "gather", nbp_glo, index_g)
1721	stream_maxresid=MAXVAL(stream_resid_glo, MASK=stream_resid_glo .LT. undef_graphfile)
1722	ENDIF
1723	CALL bcast(stream_maxresid)
1724	CALL scatter(stream_resid_glo,stream_resid_loc)
1725	stream_resid=>stream_resid_loc
1726	ENDIF
1727	!
1728	ALLOCATE (fast_reservoir(nbpt,nbasmax), stat=ier)
1729	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_reservoir','','')
1730	var_name = 'fastres'
1731	CALL ioconf_setatt_p('UNITS', 'Kg')
1732	CALL ioconf_setatt_p('LONG_NAME','Water in the fast reservoir')
1733	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g)
1734	CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero)
1735
1736	ALLOCATE (slow_reservoir(nbpt,nbasmax), stat=ier)
1737	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slow_reservoir','','')
1738	var_name = 'slowres'
1739	CALL ioconf_setatt_p('UNITS', 'Kg')
1740	CALL ioconf_setatt_p('LONG_NAME','Water in the slow reservoir')
1741	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g)
1742	CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero)
1743
1744	ALLOCATE (stream_reservoir(nbpt,nbasmax), stat=ier)
1745	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_reservoir','','')
1746	var_name = 'streamres'
1747	CALL ioconf_setatt_p('UNITS', 'Kg')
1748	CALL ioconf_setatt_p('LONG_NAME','Water in the stream reservoir')
1749	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g)
1750	CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero)
1751
1752	ALLOCATE (flood_reservoir(nbpt,nbasmax), stat=ier)
1753	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for flood_reservoir','','')
1754	var_name = 'floodres'
1755	CALL ioconf_setatt_p('UNITS', 'Kg')
1756	CALL ioconf_setatt_p('LONG_NAME','Water in the flood reservoir')
1757	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_reservoir, "gather", nbp_glo, index_g)
1758	CALL setvar_p (flood_reservoir, val_exp, 'NO_KEYWORD', zero)
1759
1760	ALLOCATE (flood_frac_bas(nbpt,nbasmax), stat=ier)
1761	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for flood_frac_bas','','')
1762	var_name = 'flood_frac_bas'
1763	CALL ioconf_setatt_p('UNITS', '-')
1764	CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per basin')
1765	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_frac_bas, "gather", nbp_glo, index_g)
1766	CALL setvar_p (flood_frac_bas, val_exp, 'NO_KEYWORD', zero)
1767
1768	ALLOCATE (flood_height(nbpt, nbasmax), stat=ier)
1769	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for flood_height','','')
1770	var_name = 'floodh'
1771	CALL ioconf_setatt_p('UNITS', '-')
1772	CALL ioconf_setatt_p('LONG_NAME','')
1773	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g)
1774	CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero)
1775
1776	ALLOCATE (pond_frac(nbpt), stat=ier)
1777	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for pond_frac','','')
1778	var_name = 'pond_frac'
1779	CALL ioconf_setatt_p('UNITS', '-')
1780	CALL ioconf_setatt_p('LONG_NAME','Pond fraction per grid box')
1781	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g)
1782	CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero)
1783
1784	var_name = 'flood_frac'
1785	CALL ioconf_setatt_p('UNITS', '-')
1786	CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per grid box')
1787	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g)
1788	CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero)
1789
1790	var_name = 'flood_res'
1791	CALL ioconf_setatt_p('UNITS','mm')
1792	CALL ioconf_setatt_p('LONG_NAME','Flooded quantity (estimation)')
1793	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_res, "gather", nbp_glo, index_g)
1794	CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero)
1795	! flood_res = zero
1796
1797	ALLOCATE (lake_reservoir(nbpt), stat=ier)
1798	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for lake_reservoir','','')
1799	var_name = 'lakeres'
1800	CALL ioconf_setatt_p('UNITS', 'Kg')
1801	CALL ioconf_setatt_p('LONG_NAME','Water in the lake reservoir')
1802	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g)
1803	CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero)
1804
1805	ALLOCATE (pond_reservoir(nbpt), stat=ier)
1806	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for pond_reservoir','','')
1807	var_name = 'pondres'
1808	CALL ioconf_setatt_p('UNITS', 'Kg')
1809	CALL ioconf_setatt_p('LONG_NAME','Water in the pond reservoir')
1810	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_reservoir, "gather", nbp_glo, index_g)
1811	CALL setvar_p (pond_reservoir, val_exp, 'NO_KEYWORD', zero)
1812	!
1813	! Map of irrigated areas
1814	!
1815	IF ( do_irrigation ) THEN
1816	ALLOCATE (irrigated(nbpt), stat=ier)
1817	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for irrigated','','')
1818	var_name = 'irrigated'
1819	CALL ioconf_setatt_p('UNITS', 'm^2')
1820	CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area')
1821	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g)
1822	CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba)
1823	ENDIF
1824
1825	IF ( do_floodplains ) THEN
1826	ALLOCATE (floodmap(nbpt), stat=ier)
1827	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodmap','','')
1828
1829	ALLOCATE (floodplains_loc(nbpt,nbasmax), stat=ier)
1830	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodplains_loc','','')
1831
1832	ALLOCATE (floodplains_glo(nbp_glo,nbasmax), stat=ier)
1833	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodplains_glo','','')
1834	IF ( .NOT. ReadGraph ) THEN
1835	IF (is_root_prc) THEN
1836	var_name = 'floodplains'
1837	CALL ioconf_setatt_p('UNITS', 'm^2')
1838	CALL ioconf_setatt_p('LONG_NAME','Surface which can be flooded')
1839	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., floodplains_glo, "gather", nbp_glo, index_g)
1840	END IF
1841	CALL scatter(floodplains_glo,floodplains_loc)
1842	floodplains=>floodplains_loc
1843	END IF
1844	ENDIF
1845	!!!
1846	!!! ANTHONY : OVERFLOW
1847	!!!
1848	IF ( dofloodoverflow) THEN
1849	ALLOCATE (orog_min_loc(nbpt,nbasmax), stat=ier)
1850	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for orog_min_loc','','')
1851	ALLOCATE (orog_min_glo(nbp_glo,nbasmax), stat=ier)
1852	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for orog_min_glo','','')
1853	!!
1854	IF ( .NOT. ReadGraph ) THEN
1855	IF (is_root_prc) THEN
1856	var_name = 'orog_min'
1857	CALL ioconf_setatt('UNITS', 'm')
1858	CALL ioconf_setatt('LONG_NAME','HTU minimum orography')
1859	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., orog_min_glo, "gather", nbp_glo, index_g)
1860	END IF
1861	CALL scatter(orog_min_glo,orog_min_loc)
1862	orog_min=>orog_min_loc
1863	END IF
1864	!
1865	ALLOCATE (route_innum_loc(nbpt,nbasmax), stat=ier)
1866	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_innum_loc','','')
1867	ALLOCATE (route_innum_glo(nbp_glo,nbasmax), stat=ier)
1868	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_innum_glo','','')
1869	CALL scatter(route_innum_glo,route_innum_loc)
1870	route_innum=>route_innum_loc
1871	!
1872	ALLOCATE (route_ingrid_loc(nbpt,nbasmax, inflows), stat=ier)
1873	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_ingrid_loc','','')
1874	ALLOCATE (route_ingrid_glo(nbp_glo,nbasmax,inflows), stat=ier)
1875	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_ingrid_glo','','')
1876	CALL scatter(route_ingrid_glo,route_ingrid_loc)
1877	route_ingrid=>route_ingrid_loc
1878	!
1879	ALLOCATE (route_inbasin_loc(nbpt,nbasmax, inflows), stat=ier)
1880	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_inbasin_loc','','')
1881	ALLOCATE (route_inbasin_glo(nbp_glo,nbasmax, inflows), stat=ier)
1882	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_inbasin_glo','','')
1883	CALL scatter(route_inbasin_glo,route_inbasin_loc)
1884	route_inbasin=>route_inbasin_loc
1885	END IF
1886	!!!
1887	!!!
1888	!!!
1889	IF ( doswamps ) THEN
1890	ALLOCATE (swamp(nbpt), stat=ier)
1891	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for swamp','','')
1892	var_name = 'swamp'
1893	CALL ioconf_setatt_p('UNITS', 'm^2')
1894	CALL ioconf_setatt_p('LONG_NAME','Surface which can become swamp')
1895	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., swamp, "gather", nbp_glo, index_g)
1896	CALL setvar_p (swamp, val_exp, 'NO_KEYWORD', undef_sechiba)
1897	ENDIF
1898	!
1899	! Put into the restart file the fluxes so that they can be regenerated at restart.
1900	!
1901	ALLOCATE (lakeinflow_mean(nbpt), stat=ier)
1902	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for lakeinflow_mean','','')
1903	var_name = 'lakeinflow'
1904	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
1905	CALL ioconf_setatt_p('LONG_NAME','Lake inflow')
1906	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
1907	CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
1908
1909	ALLOCATE (returnflow_mean(nbpt), stat=ier)
1910	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for returnflow_mean','','')
1911	var_name = 'returnflow'
1912	CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
1913	CALL ioconf_setatt_p('LONG_NAME','Deep return flux')
1914	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g)
1915	CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero)
1916	returnflow(:) = returnflow_mean(:)
1917
1918	ALLOCATE (reinfiltration_mean(nbpt), stat=ier)
1919	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for reinfiltration_mean','','')
1920	var_name = 'reinfiltration'
1921	CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
1922	CALL ioconf_setatt_p('LONG_NAME','Top return flux')
1923	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinfiltration_mean, "gather", nbp_glo, index_g)
1924	CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero)
1925	reinfiltration(:) = reinfiltration_mean(:)
1926
1927	ALLOCATE (irrigation_mean(nbpt), stat=ier)
1928	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for irrigation_mean','','')
1929	ALLOCATE (irrig_netereq(nbpt), stat=ier)
1930	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for irrig_netereq','','')
1931	irrig_netereq(:) = zero
1932
1933	IF ( do_irrigation ) THEN
1934	var_name = 'irrigation'
1935	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
1936	CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux')
1937	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g)
1938	CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero)
1939	ELSE
1940	irrigation_mean(:) = zero
1941	ENDIF
1942	irrigation(:) = irrigation_mean(:)
1943
1944	ALLOCATE (fast_temp(nbpt,nbasmax), stat=ier)
1945	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_temp','','')
1946	var_name = 'fasttemp'
1947	CALL ioconf_setatt_p('UNITS', 'K')
1948	CALL ioconf_setatt_p('LONG_NAME','Water temperature in the fast reservoir')
1949	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_temp, "gather", nbp_glo, index_g)
1950
1951	ALLOCATE (slow_temp(nbpt,nbasmax), stat=ier)
1952	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slow_temp','','')
1953	var_name = 'slowtemp'
1954	CALL ioconf_setatt_p('UNITS', 'K')
1955	CALL ioconf_setatt_p('LONG_NAME','Water temperature in the slow reservoir')
1956	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_temp, "gather", nbp_glo, index_g)
1957
1958	IF ( COUNT(fast_temp == val_exp) == nbpt*nbasmax ) THEN
1959	CALL groundwatertemp(nbpt, nbasmax, ngrnd, tempdiag, znt, dlt, fast_temp, slow_temp)
1960	ENDIF
1961
1962	ALLOCATE (stream_temp(nbpt,nbasmax), stat=ier)
1963	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_temp','','')
1964	var_name = 'streamtemp'
1965	CALL ioconf_setatt_p('UNITS', 'K')
1966	CALL ioconf_setatt_p('LONG_NAME','Water temperature in the stream reservoir')
1967	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_temp, "gather", nbp_glo, index_g)
1968
1969	IF ( COUNT(stream_temp == val_exp) == nbpt*nbasmax ) THEN
1970	DO ig=1,nbpt
1971	stream_temp(ig,:) = tempdiag(ig,1)
1972	ENDDO
1973	ENDIF
1974
1975	ALLOCATE (riverflow_mean(nbpt), stat=ier)
1976	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for riverflow_mean','','')
1977	var_name = 'riverflow'
1978	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
1979	CALL ioconf_setatt_p('LONG_NAME','River flux into the sea')
1980	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
1981	CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
1982	riverflow(:) = riverflow_mean(:)
1983
1984	ALLOCATE (coastalflow_mean(nbpt), stat=ier)
1985	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for coastalflow_mean','','')
1986	var_name = 'coastalflow'
1987	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
1988	CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea')
1989	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
1990	CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
1991	coastalflow(:) = coastalflow_mean(:)
1992
1993	! Locate it at the 2m level
1994	ipn = MINLOC(ABS(zlt-2))
1995	floodtemp_lev = ipn(1)
1996	ALLOCATE (floodtemp(nbpt), stat=ier)
1997	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodtemp','','')
1998	floodtemp(:) = tempdiag(:,floodtemp_lev)
1999
2000	ALLOCATE(hydrographs(nbpt), stat=ier)
2001	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for hydrographs','','')
2002	var_name = 'hydrographs'
2003	CALL ioconf_setatt_p('UNITS', 'kg/dt_sechiba')
2004	CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid')
2005	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g)
2006	CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero)
2007
2008	ALLOCATE(hydrotemp(nbpt), stat=ier)
2009	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for hydrotemp','','')
2010	var_name = 'hydrotemp'
2011	CALL ioconf_setatt_p('UNITS', 'K')
2012	CALL ioconf_setatt_p('LONG_NAME','Temperature of most significant river of grid')
2013	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrotemp, "gather", nbp_glo, index_g)
2014	CALL setvar_p (hydrotemp, val_exp, 'NO_KEYWORD', ZeroCelsius)
2015
2016	ALLOCATE(slowflow_diag(nbpt), stat=ier)
2017	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slowflow_diag','','')
2018	var_name = 'slowflow_diag'
2019	CALL ioconf_setatt_p('UNITS', 'kg/dt_sechiba')
2020	CALL ioconf_setatt_p('LONG_NAME','Slowflow hydrograph at outlow of grid')
2021	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,slowflow_diag, "gather", nbp_glo, index_g)
2022	CALL setvar_p (slowflow_diag, val_exp, 'NO_KEYWORD', zero)
2023	!
2024	! Grid diagnostic at representative HTU
2025	!
2026	ALLOCATE(hydrodiag_loc(nbpt),hydrodiag_glo(nbp_glo),stat=ier)
2027	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for hydrodiag_glo','','')
2028	IF ( .NOT. ReadGraph ) THEN
2029	IF (is_root_prc) THEN
2030	ALLOCATE(tmp_real(nbp_glo))
2031	var_name = 'gridrephtu'
2032	CALL ioconf_setatt('UNITS', '-')
2033	CALL ioconf_setatt('LONG_NAME','Representative HTU for the grid')
2034	CALL restget(rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,tmp_real, "gather", nbp_glo, index_g)
2035	hydrodiag_glo(:) = 1
2036	WHERE ( tmp_real .LT. val_exp )
2037	hydrodiag_glo = NINT(tmp_real)
2038	ENDWHERE
2039	DEALLOCATE(tmp_real)
2040	ENDIF
2041	CALL scatter(hydrodiag_glo, hydrodiag_loc)
2042	ENDIF
2043	!
2044	! Station diagnostics
2045	!
2046	ALLOCATE(HTUdiag_loc(nbpt,nbasmon), stat=ier)
2047	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUdiag_loc','','')
2048	ALLOCATE(HTUdiag_glo(nbp_glo,nbasmon), stat=ier)
2049	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUdiag_glo','','')
2050	ALLOCATE(tmp_real_g2(nbp_glo,nbasmon), stat=ier)
2051	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tmp_real_g2','','')
2052	!
2053	IF (is_root_prc) THEN
2054	var_name = 'htudiag'
2055	CALL ioconf_setatt('UNITS', 'index')
2056	CALL ioconf_setatt('LONG_NAME','Index of HTU to be monitored')
2057	CALL restget(rest_id, var_name, nbp_glo, nbasmon, 1, kjit, .TRUE., tmp_real_g2, "gather", nbp_glo, index_g)
2058	HTUdiag_glo(:,:) = -1
2059	WHERE ( tmp_real_g2 .LT. val_exp )
2060	HTUdiag_glo = NINT(tmp_real_g2)
2061	ENDWHERE
2062	nbhtumon = 0
2063	DO ig=1,nbp_glo
2064	DO im=1,nbasmon
2065	IF ( HTUdiag_glo(ig,im) > 0 ) THEN
2066	nbhtumon = nbhtumon + 1
2067	ENDIF
2068	ENDDO
2069	ENDDO
2070	WRITE(numout,*) "After restget : Found a total of ", nbhtumon, " HTUs to be monitored and written into HTUhgmon"
2071	ENDIF
2072	CALL scatter(HTUdiag_glo, HTUdiag_loc)
2073	CALL bcast(nbhtumon)
2074	DEALLOCATE(tmp_real_g2)
2075	!
2076	ALLOCATE(HTUhgmon(nbpt,nbasmon), stat=ier)
2077	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUhgmon','','')
2078	ALLOCATE(HTUhgmon_glo(nbp_glo,nbasmon), stat=ier)
2079	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUhgmon_glo','','')
2080	!
2081	IF (is_root_prc) THEN
2082	var_name = 'htuhgmon'
2083	CALL ioconf_setatt('UNITS', 'kg/dt_sechiba')
2084	CALL ioconf_setatt('LONG_NAME','Hydrograph at selected HTU of grid')
2085	CALL restget(rest_id, var_name, nbp_glo, nbasmon, 1, kjit, .TRUE., HTUhgmon_glo, "gather", nbp_glo, index_g)
2086	WHERE ( HTUhgmon_glo .GE. val_exp )
2087	HTUhgmon_glo = zero
2088	ENDWHERE
2089	ENDIF
2090	CALL scatter(HTUhgmon_glo, HTUhgmon)
2091	DEALLOCATE(HTUhgmon_glo)
2092	!
2093	! Restart of the temperature monitoring
2094	!
2095	ALLOCATE(HTUtempmon(nbpt,nbasmon), stat=ier)
2096	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUtempmon','','')
2097	ALLOCATE(HTUtempmon_glo(nbp_glo,nbasmon), stat=ier)
2098	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUtempmon_glo','','')
2099	!
2100	IF (is_root_prc) THEN
2101	var_name = 'htutempmon'
2102	CALL ioconf_setatt('UNITS', 'K')
2103	CALL ioconf_setatt('LONG_NAME','Temperature at selected HTU of grid')
2104	CALL restget(rest_id, var_name, nbp_glo, nbasmon, 1, kjit, .TRUE., HTUtempmon_glo, "gather", nbp_glo, index_g)
2105	WHERE ( HTUtempmon_glo .GE. val_exp )
2106	HTUtempmon_glo = ZeroCelsius
2107	ENDWHERE
2108	HTUtempmon_glo(:,:) = ZeroCelsius
2109	ENDIF
2110	CALL scatter(HTUtempmon_glo, HTUtempmon)
2111	DEALLOCATE(HTUtempmon_glo)
2112	!
2113	! The diagnostic variables, they are initialized from the above restart variables.
2114	!
2115	ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), flood_diag(nbpt), &
2116	& pond_diag(nbpt), lake_diag(nbpt), stat=ier)
2117	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_diag,..','','')
2118
2119	fast_diag(:) = zero
2120	slow_diag(:) = zero
2121	stream_diag(:) = zero
2122	flood_diag(:) = zero
2123	pond_diag(:) = zero
2124	lake_diag(:) = zero
2125
2126	DO ig=1,nbpt
2127	totarea = zero
2128	DO ib=1,nbasmax
2129	totarea = totarea + routing_area(ig,ib)
2130	fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib)
2131	slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib)
2132	stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib)
2133	flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib)
2134	ENDDO
2135	!
2136	fast_diag(ig) = fast_diag(ig)/totarea
2137	slow_diag(ig) = slow_diag(ig)/totarea
2138	stream_diag(ig) = stream_diag(ig)/totarea
2139	flood_diag(ig) = flood_diag(ig)/totarea
2140	!
2141	! This is the volume of the lake scaled to the entire grid.
2142	! It would be better to scale it to the size of the lake
2143	! but this information is not yet available.
2144	!
2145	lake_diag(ig) = lake_reservoir(ig)/totarea
2146	!
2147	ENDDO
2148	!
2149	! Get from the restart the fluxes we accumulated.
2150	!
2151	ALLOCATE (floodout_mean(nbpt), stat=ier)
2152	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodout_mean','','')
2153	var_name = 'floodout_route'
2154	CALL ioconf_setatt_p('UNITS', 'Kg')
2155	CALL ioconf_setatt_p('LONG_NAME','Accumulated flow out of floodplains for routing')
2156	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g)
2157	CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero)
2158
2159	ALLOCATE (runoff_mean(nbpt), stat=ier)
2160	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for runoff_mean','','')
2161	var_name = 'runoff_route'
2162	CALL ioconf_setatt_p('UNITS', 'Kg')
2163	CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing')
2164	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
2165	CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
2166
2167	ALLOCATE(drainage_mean(nbpt), stat=ier)
2168	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for drainage_mean','','')
2169	var_name = 'drainage_route'
2170	CALL ioconf_setatt_p('UNITS', 'Kg')
2171	CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing')
2172	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
2173	CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
2174
2175	ALLOCATE(transpot_mean(nbpt), stat=ier)
2176	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for transpot_mean','','')
2177	var_name = 'transpot_route'
2178	CALL ioconf_setatt_p('UNITS', 'Kg/m^2')
2179	CALL ioconf_setatt_p('LONG_NAME','Accumulated potential transpiration for routing/irrigation')
2180	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., transpot_mean, "gather", nbp_glo, index_g)
2181	CALL setvar_p (transpot_mean, val_exp, 'NO_KEYWORD', zero)
2182
2183	ALLOCATE(precip_mean(nbpt), stat=ier)
2184	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for precip_mean','','')
2185	var_name = 'precip_route'
2186	CALL ioconf_setatt_p('UNITS', 'Kg/m^2')
2187	CALL ioconf_setatt_p('LONG_NAME','Accumulated rain precipitation for irrigation')
2188	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g)
2189	CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero)
2190
2191	ALLOCATE(humrel_mean(nbpt), stat=ier)
2192	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for humrel_mean','','')
2193	var_name = 'humrel_route'
2194	CALL ioconf_setatt_p('UNITS', '-')
2195	CALL ioconf_setatt_p('LONG_NAME','Mean humrel for irrigation')
2196	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g)
2197	CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un)
2198
2199	ALLOCATE(k_litt_mean(nbpt), stat=ier)
2200	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for k_litt_mean','','')
2201	var_name = 'k_litt_route'
2202	CALL ioconf_setatt_p('UNITS', '-')
2203	CALL ioconf_setatt_p('LONG_NAME','Mean cond. for litter')
2204	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g)
2205	CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero)
2206
2207	ALLOCATE(totnobio_mean(nbpt), stat=ier)
2208	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for totnobio_mean','','')
2209	var_name = 'totnobio_route'
2210	CALL ioconf_setatt_p('UNITS', '-')
2211	CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of no bio for irrigation')
2212	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g)
2213	CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero)
2214
2215	ALLOCATE(vegtot_mean(nbpt), stat=ier)
2216	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for vegtot_mean','','')
2217	var_name = 'vegtot_route'
2218	CALL ioconf_setatt_p('UNITS', '-')
2219	CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of vegetation')
2220	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g)
2221	CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un)
2222	!
2223	ALLOCATE(tempdiag_mean(nbpt,ngrnd), stat=ier)
2224	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tempdiag_mean','','')
2225	var_name = 'tempdiag_route'
2226	CALL ioconf_setatt_p('UNITS', 'K')
2227	CALL ioconf_setatt_p('LONG_NAME','Mean temperature profile')
2228	CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., tempdiag_mean, "gather", nbp_glo, index_g)
2229	CALL setvar_p (tempdiag_mean, val_exp, 'NO_KEYWORD', Zero)
2230	!
2231	DEALLOCATE(tmp_real_g)
2232	!
2233	! Other variables
2234	!
2235	ALLOCATE(streamlimit(nbpt), stat=ier)
2236	!
2237	END SUBROUTINE routing_hr_init
2238	!
2239	!! ================================================================================================================================
2240	!! SUBROUTINE : routing_hr_restartconsistency
2241	!!
2242	!>\BRIEF : This subroutine will verify that the important dimensions for the routing exist in the restart file
2243	!! or can be read from the routing_graph.nc file. This ensures that if the restart and graph file are
2244	!! not consistent the latter is read and overwrite whatever was in the restart. Then the user will be
2245	!! using the new routing graph from the routing_graph.nc and not whatever is in the restart.
2246	!! \n
2247	!_ ================================================================================================================================
2248
2249	SUBROUTINE routing_hr_restartconsistency(varname, dimgraph, dimrestart)
2250	!
2251	IMPLICIT NONE
2252	!
2253	!! INPUT VARIABLES
2254	CHARACTER(LEN=80), INTENT(in) :: varname !! Name of dimension
2255	INTEGER(i_std), INTENT(inout) :: dimgraph !! Dimension read in the graph file and also final result.
2256	REAL(r_std), INTENT(in) :: dimrestart !! Dimension read in the restart file.
2257	!
2258	! Case when the routing_hr_graphinfo could not get any information from the routing_graph.nc
2259	IF ( dimgraph < un ) THEN
2260	IF ( dimrestart > zero ) THEN
2261	! Only information from the restart.
2262	dimgraph = NINT(dimrestart)
2263	ELSE
2264	WRITE(,) "Problem : No information in the routing_graph file and no routing information in restart ", TRIM(varname)
2265	CALL ipslerr(3,'routing_hr_restartconsistency',&
2266	'No routing_graph file availble and no information in restart.', &
2267	'Cannot perform routing in ORCHIDEE.', ' ')
2268	ENDIF
2269
2270	! Information from the routing_graph.nc file exists !
2271	ELSE
2272	IF ( dimgraph .NE. NINT(dimrestart) ) THEN
2273	WRITE(,) "Problem for ", TRIM(varname)," in restart is not the same as in routing_graph.nc "
2274	WRITE(,) "Value of ", TRIM(varname), " in restart file : ", dimrestart
2275	WRITE(,) "Value of ", TRIM(varname), " in routing_graph.nc file : ", dimgraph
2276	CALL ipslerr(2,'routing_hr_restartconsistency',&
2277	'The value of dimension provided is not consistant with the one in routing_graph file.', &
2278	'We will read a new graph from the given file.', ' ')
2279	ReadGraph = .TRUE.
2280	ELSE
2281	!! Nothing to do
2282	ENDIF
2283	ENDIF
2284	END SUBROUTINE routing_hr_restartconsistency
2285	!! ================================================================================================================================
2286	!! SUBROUTINE : routing_highres_clear
2287	!!
2288	!>\BRIEF : This subroutine deallocates the block memory previously allocated.
2289	!! \n
2290	!_ ================================================================================================================================
2291
2292	SUBROUTINE routing_highres_clear()
2293
2294	IF (ALLOCATED(routing_area_loc)) DEALLOCATE(routing_area_loc)
2295	IF (ALLOCATED(route_togrid_loc)) DEALLOCATE(route_togrid_loc)
2296	IF (ALLOCATED(route_tobasin_loc)) DEALLOCATE(route_tobasin_loc)
2297	IF (ALLOCATED(route_nbintobas_loc)) DEALLOCATE(route_nbintobas_loc)
2298	IF (ALLOCATED(global_basinid_loc)) DEALLOCATE(global_basinid_loc)
2299	IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc)
2300	IF (ALLOCATED(stream_resid_loc)) DEALLOCATE(stream_resid_loc)
2301	IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo)
2302	IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo)
2303	IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo)
2304	IF (ALLOCATED(route_nbintobas_glo)) DEALLOCATE(route_nbintobas_glo)
2305	IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo)
2306	IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo)
2307	IF (ALLOCATED(stream_resid_glo)) DEALLOCATE(stream_resid_glo)
2308	IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir)
2309	IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir)
2310	IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir)
2311
2312	IF (ALLOCATED(fast_temp)) DEALLOCATE(fast_temp)
2313	IF (ALLOCATED(slow_temp)) DEALLOCATE(slow_temp)
2314	IF (ALLOCATED(stream_temp)) DEALLOCATE(stream_temp)
2315
2316	IF (ALLOCATED(flood_reservoir)) DEALLOCATE(flood_reservoir)
2317	IF (ALLOCATED(flood_frac_bas)) DEALLOCATE(flood_frac_bas)
2318	IF (ALLOCATED(flood_height)) DEALLOCATE(flood_height)
2319	IF (ALLOCATED(pond_frac)) DEALLOCATE(pond_frac)
2320	IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir)
2321	IF (ALLOCATED(pond_reservoir)) DEALLOCATE(pond_reservoir)
2322	IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_mean)
2323	IF (ALLOCATED(reinfiltration_mean)) DEALLOCATE(reinfiltration_mean)
2324	IF (ALLOCATED(riverflow_mean)) DEALLOCATE(riverflow_mean)
2325	IF (ALLOCATED(coastalflow_mean)) DEALLOCATE(coastalflow_mean)
2326	IF (ALLOCATED(lakeinflow_mean)) DEALLOCATE(lakeinflow_mean)
2327	IF (ALLOCATED(runoff_mean)) DEALLOCATE(runoff_mean)
2328	IF (ALLOCATED(floodout_mean)) DEALLOCATE(floodout_mean)
2329	IF (ALLOCATED(drainage_mean)) DEALLOCATE(drainage_mean)
2330	IF (ALLOCATED(transpot_mean)) DEALLOCATE(transpot_mean)
2331	IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean)
2332	IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean)
2333	IF (ALLOCATED(k_litt_mean)) DEALLOCATE(k_litt_mean)
2334	IF (ALLOCATED(tempdiag_mean)) DEALLOCATE(tempdiag_mean)
2335	IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean)
2336	IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean)
2337	IF (ALLOCATED(floodtemp)) DEALLOCATE(floodtemp)
2338	IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc)
2339	IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo)
2340	IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs)
2341	IF (ALLOCATED(hydrotemp)) DEALLOCATE(hydrotemp)
2342	IF (ALLOCATED(HTUhgmon)) DEALLOCATE(HTUhgmon)
2343	IF (ALLOCATED(HTUhgmon_glo)) DEALLOCATE(HTUhgmon_glo)
2344	IF (ALLOCATED(HTUtempmon)) DEALLOCATE(HTUtempmon)
2345	IF (ALLOCATED(HTUtempmon_glo)) DEALLOCATE(HTUtempmon_glo)
2346	IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag)
2347	IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean)
2348	IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated)
2349	IF (ALLOCATED(floodplains_glo)) DEALLOCATE(floodplains_glo)
2350	IF (ALLOCATED(floodplains_loc)) DEALLOCATE(floodplains_loc)
2351	IF (ALLOCATED(swamp)) DEALLOCATE(swamp)
2352	IF (ALLOCATED(fast_diag)) DEALLOCATE(fast_diag)
2353	IF (ALLOCATED(slow_diag)) DEALLOCATE(slow_diag)
2354	IF (ALLOCATED(stream_diag)) DEALLOCATE(stream_diag)
2355	IF (ALLOCATED(flood_diag)) DEALLOCATE(flood_diag)
2356	IF (ALLOCATED(pond_diag)) DEALLOCATE(pond_diag)
2357	IF (ALLOCATED(lake_diag)) DEALLOCATE(lake_diag)
2358	!
2359	IF (ALLOCATED(route_innum_loc)) DEALLOCATE(route_innum_loc)
2360	IF (ALLOCATED(route_ingrid_loc)) DEALLOCATE(route_ingrid_loc)
2361	IF (ALLOCATED(route_inbasin_loc)) DEALLOCATE(route_inbasin_loc)
2362	IF (ALLOCATED(route_innum_glo)) DEALLOCATE(route_innum_glo)
2363	IF (ALLOCATED(route_ingrid_glo)) DEALLOCATE(route_ingrid_glo)
2364	IF (ALLOCATED(route_inbasin_glo)) DEALLOCATE(route_inbasin_glo)
2365	!
2366	IF (ALLOCATED(orog_min_loc)) DEALLOCATE(orog_min_loc)
2367	IF (ALLOCATED(orog_min_glo)) DEALLOCATE(orog_min_glo)
2368	!
2369	IF (ALLOCATED(floodcri_loc)) DEALLOCATE(floodcri_loc)
2370	IF (ALLOCATED(floodcri_glo)) DEALLOCATE(floodcri_glo)
2371	IF (ALLOCATED(fp_beta_loc)) DEALLOCATE(fp_beta_loc)
2372	IF (ALLOCATED(fp_beta_glo)) DEALLOCATE(fp_beta_glo)
2373
2374	END SUBROUTINE routing_highres_clear
2375	!
2376
2377	!! ================================================================================================================================
2378	!! SUBROUTINE : routing_hr_flow
2379	!!
2380	!>\BRIEF This subroutine computes the transport of water in the various reservoirs
2381	!! (including ponds and floodplains) and the water withdrawals from the reservoirs for irrigation.
2382	!!
2383	!! DESCRIPTION (definitions, functional, design, flags) :
2384	!! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an
2385	!! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes
2386	!! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of
2387	!! the up-stream basin and adding it to the down-stream one.
2388	!! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once
2389	!! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them
2390	!! the hydrographs of the largest rivers we have chosen to monitor.
2391	!!
2392	!! RECENT CHANGE(S): None
2393	!!
2394	!! MAIN OUTPUT VARIABLE(S): lakeinflow, returnflow, reinfiltration, irrigation, riverflow, coastalflow, hydrographs, flood_frac, flood_res
2395	!!
2396	!! REFERENCES :
2397	!! - Ngo-Duc, T., K. Laval, G. Ramillien, J. Polcher, and A. Cazenave (2007)
2398	!! Validation of the land water storage simulated by Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) with Gravity Recovery and Climate Experiment (GRACE) data.
2399	!! Water Resour. Res., 43, W04427, doi:10.1029/2006WR004941.
2400	!! * Irrigation:
2401	!! - de Rosnay, P., J. Polcher, K. Laval, and M. Sabre (2003)
2402	!! Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula.
2403	!! Geophys. Res. Lett., 30(19), 1986, doi:10.1029/2003GL018024.
2404	!! - A.C. Vivant (2003)
2405	!! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte.
2406	!! , , 51pp.
2407	!! - N. Culson (2004)
2408	!! Impact de l'irrigation sur le cycle de l'eau
2409	!! Master thesis, Paris VI University, 55pp.
2410	!! - X.-T. Nguyen-Vinh (2005)
2411	!! Analyse de l'impact de l'irrigation en Amerique du Nord - plaine du Mississippi - sur la climatologie regionale
2412	!! Master thesis, Paris VI University, 33pp.
2413	!! - M. Guimberteau (2006)
2414	!! Analyse et modifications proposees de la modelisation de l'irrigation dans un modele de surface.
2415	!! Master thesis, Paris VI University, 46pp.
2416	!! - Guimberteau M. (2010)
2417	!! Modelisation de l'hydrologie continentale et influences de l'irrigation sur le cycle de l'eau.
2418	!! Ph.D. thesis, Paris VI University, 195pp.
2419	!! - Guimberteau M., Laval K., Perrier A. and Polcher J. (2011).
2420	!! Global effect of irrigation and its impact on the onset of the Indian summer monsoon.
2421	!! In press, Climate Dynamics, doi: 10.1007/s00382-011-1252-5.
2422	!! * Floodplains:
2423	!! - A.C. Vivant (2002)
2424	!! L'ecoulement lateral de l'eau sur les surfaces continentales. Prise en compte des plaines d'inondations dans ORCHIDEE.
2425	!! Master thesis, Paris VI University, 46pp.
2426	!! - A.C. Vivant (2003)
2427	!! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte.
2428	!! , , 51pp.
2429	!! - T. d'Orgeval (2006)
2430	!! Impact du changement climatique sur le cycle de l'eau en Afrique de l'Ouest: modelisation et incertitudes.
2431	!! Ph.D. thesis, Paris VI University, 188pp.
2432	!! - T. d'Orgeval, J. Polcher, and P. de Rosnay (2008)
2433	!! Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes.
2434	!! Hydrol. Earth Syst. Sci., 12, 1387-1401
2435	!! - M. Guimberteau, G. Drapeau, J. Ronchail, B. Sultan, J. Polcher, J.-M. Martinez, C. Prigent, J.-L. Guyot, G. Cochonneau,
2436	!! J. C. Espinoza, N. Filizola, P. Fraizy, W. Lavado, E. De Oliveira, R. Pombosa, L. Noriega, and P. Vauchel (2011)
2437	!! Discharge simulation in the sub-basins of the Amazon using ORCHIDEE forced by new datasets.
2438	!! Hydrol. Earth Syst. Sci. Discuss., 8, 11171-11232, doi:10.5194/hessd-8-11171-2011
2439	!!
2440	!! FLOWCHART :None
2441	!! \n
2442	!_ ================================================================================================================================
2443
2444	SUBROUTINE routing_hr_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, &
2445	& vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, tempdiag, &
2446	& reinf_slope, lakeinflow, returnflow, reinfiltration, irrigation, riverflow, &
2447	& coastalflow, hydrographs, slowflow_diag, flood_frac, flood_res, &
2448	& netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, &
2449	& stemp_total_tend, stemp_advec_tend, stemp_relax_tend)
2450	!
2451	IMPLICIT NONE
2452	!
2453	!! INPUT VARIABLES
2454	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
2455	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
2456	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
2457	REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff (kg/m^2/dt)
2458	REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt)
2459	REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt)
2460	REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated fraction (unitless;0-1)
2461	REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas which can not have vegetation
2462	REAL(r_std), INTENT(in) :: transpot_mean(nbpt) !! Mean potential transpiration of the vegetation (kg/m^2/dt)
2463	REAL(r_std), INTENT(in) :: precip(nbpt) !! Rainfall (kg/m^2/dt)
2464	REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless)
2465	REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
2466	REAL(r_std), INTENT(in) :: floodtemp(nbpt) !! Temperature to decide if floodplains work (K)
2467	REAL(r_std), INTENT(in) :: tempdiag(nbpt,ngrnd) !! Soil temperature profiles (K)
2468	REAL(r_std), INTENT(in) :: reinf_slope(nbpt) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
2469	REAL(r_std), INTENT(out) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt)
2470	!
2471	!! OUTPUT VARIABLES
2472	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns into the grid box.
2473	!! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt_routing)
2474	REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
2475	REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt_routing)
2476	REAL(r_std), INTENT(out) :: riverflow(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt_routing)
2477	REAL(r_std), INTENT(out) :: coastalflow(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt_routing)
2478	REAL(r_std), INTENT(out) :: hydrographs(nbpt) !! Hydrographs at the outflow of the grid box for major basins (kg/dt)
2479	REAL(r_std), INTENT(out) :: slowflow_diag(nbpt) !! Hydrographs of slow_flow = routed slow_flow for major basins (kg/dt)
2480	REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
2481	REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
2482
2483	REAL(r_std), INTENT(out) :: netflow_stream_diag(nbpt) !! Input - Output flow to stream reservoir
2484	REAL(r_std), INTENT(out) :: netflow_fast_diag(nbpt) !! Input - Output flow to fast reservoir
2485	REAL(r_std), INTENT(out) :: netflow_slow_diag(nbpt) !! Input - Output flow to slow reservoir
2486	REAL(r_std), INTENT(out) :: stemp_total_tend(nbpt, nbasmax) !! Total tendency in GJ/s computed for the stream reservoir.
2487	REAL(r_std), INTENT(out) :: stemp_advec_tend(nbpt, nbasmax) !! Tendency (GJ/s) produced by advection
2488	REAL(r_std), INTENT(out) :: stemp_relax_tend(nbpt, nbasmax) !! Tendency (GJ/s) produced by relaxation
2489	!
2490	!! LOCAL VARIABLES
2491	REAL(r_std), DIMENSION(nbpt, nbasmax) :: fast_flow !! Outflow from the fast reservoir (kg/dt)
2492	REAL(r_std), DIMENSION(nbpt, nbasmax) :: slow_flow !! Outflow from the slow reservoir (kg/dt)
2493	REAL(r_std), DIMENSION(nbpt, nbasmax) :: stream_flow !! Outflow from the stream reservoir (kg/dt)
2494	REAL(r_std), DIMENSION(nbpt, nbasmax) :: flood_flow !! Outflow from the floodplain reservoir (kg/dt)
2495	REAL(r_std), DIMENSION(nbpt, nbasmax) :: pond_inflow !! Inflow to the pond reservoir (kg/dt)
2496	REAL(r_std), DIMENSION(nbpt, nbasmax) :: pond_drainage !! Drainage from pond (kg/m^2/dt)
2497	REAL(r_std), DIMENSION(nbpt, nbasmax) :: flood_drainage !! Drainage from floodplains (kg/m^2/dt)
2498	REAL(r_std), DIMENSION(nbpt, nbasmax) :: return_swamp !! Inflow to the swamp (kg/dt)
2499	REAL(r_std), DIMENSION(nbpt, nbasmax) :: source
2500	REAL(r_std), DIMENSION(nbpt, nbasmax) :: ewh
2501	!
2502	! Irrigation per basin
2503	!
2504	REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_needs !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
2505	REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_actual !! Possible irrigation according to the water availability in the reservoirs (kg)
2506	REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_deficit !! Amount of water missing for irrigation (kg)
2507	REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_adduct !! Amount of water carried over from other basins for irrigation (kg)
2508	!
2509	! The transport terms are over a larger indexing space so that outlfows to ocean and lakes do not generate out of bounds issues.
2510	! Non existing HTU have their index set to zero and their memory will end-up in index 0 of transport.
2511	!
2512	REAL(r_std), DIMENSION(nbpt, 0:nbasmax+3) :: transport !! Water transport between basins (kg/dt)
2513	REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_glo !! Water transport between basins (kg/dt)
2514	REAL(r_std), DIMENSION(nbpt, 0:nbasmax+3) :: transport_temp !! Temperature transport between grids
2515	REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_temp_glo !! Temperature transport global for transfers
2516	!
2517	REAL(r_std) :: oldtemp
2518	REAL(r_std) :: oldstream
2519	INTEGER(i_std), SAVE :: nbunpy=0
2520	!
2521	REAL(r_std), DIMENSION(nbpt, nbasmax) :: floods !! Water flow in to the floodplains (kg/dt)
2522	REAL(r_std), DIMENSION(nbpt, nbasmax) :: potflood !! Potential inflow to the swamps (kg/dt)
2523	REAL(r_std), DIMENSION(nbpt) :: tobeflooded !! Maximal surface which can be inundated in each grid box (m^2)
2524	REAL(r_std), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2)
2525	REAL(r_std), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg)
2526	REAL(r_std), DIMENSION(nbasmax) :: pond_excessflow !!
2527	REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt)
2528	REAL(r_std) :: floodindex !! Fraction of grid box area inundated (unitless;0-1)
2529	REAL(r_std) :: pondex !!
2530	REAL(r_std) :: stream_tot !! Total water amount in the stream reservoirs (kg)
2531	REAL(r_std) :: adduction !! Importation of water from a stream reservoir of a neighboring grid box (kg)
2532	REAL(r_std), DIMENSION(nbp_glo) :: lake_overflow_g !! Removed water from lake reservoir on global grid (kg/gridcell/dt_routing)
2533	REAL(r_std), DIMENSION(nbpt) :: lake_overflow !! Removed water from lake reservoir on local grid (kg/gridcell/dt_routing)
2534	REAL(r_std), DIMENSION(nbpt) :: lake_overflow_coast !! lake_overflow distributed on coast gridcells, only diag(kg/gridcell/dt_routing)
2535	REAL(r_std) :: total_lake_overflow !! Sum of lake_overflow over full grid (kg)
2536	REAL(r_std), DIMENSION(8,nbasmax) :: streams_around !! Stream reservoirs of the neighboring grid boxes (kg)
2537	INTEGER(i_std), DIMENSION(8) :: igrd !!
2538	INTEGER(i_std), DIMENSION(2) :: ff !!
2539	INTEGER(i_std), DIMENSION(1) :: fi !!
2540	INTEGER(i_std) :: ig, ib, ib2, ig2, im !! Indices (unitless)
2541	INTEGER(i_std) :: rtg, rtb, in, ing, inb,inf!! Indices (unitless)
2542	!INTEGER(i_std) :: numflood !!
2543	INTEGER(i_std) :: ier, negslow !! Error handling
2544	INTEGER(i_std), DIMENSION(20) :: negig, negib
2545	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_flow_g !! Outflow from the fast reservoir (kg/dt)
2546	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: slow_flow_g !! Outflow from the slow reservoir (kg/dt)
2547	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stream_flow_g !! Outflow from the stream reservoir (kg/dt)
2548	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_temp_g !! Temperature of the fast reservoir (K)
2549	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: slow_temp_g !! Temperature of the slow reservoir (K)
2550	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stream_temp_g !! Temperature of the stream reservoir (K)
2551	!REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_height_g !! Floodplains height (m)
2552	!REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas_g !! Fraction of the HTU flooded
2553	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrig_deficit_glo !! Amount of water missing for irrigation (kg)
2554	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stream_reservoir_glo !! Water amount in the stream reservoir (kg)
2555	!REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir_glo !! Water amount in the stream reservoir (kg)
2556	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrig_adduct_glo !! Amount of water carried over from other basins for irrigation (kg)
2557
2558	REAL(r_std) :: reduced !! Discharge reduction due to floodplains
2559	REAL(r_std) :: htmp, hscale !! Water height scalingfor temperature relaxation
2560	REAL(r_std) :: krelax, den
2561	!! PARAMETERS
2562	LOGICAL, PARAMETER :: check_reservoir = .FALSE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false)
2563	!_ ================================================================================================================================
2564	!
2565	!
2566	hscale = 1.
2567	CALL getin_p('ROUTING_HSCALEKH',hscale)
2568	!
2569	transport(:,:) = zero
2570	transport_glo(:,:) = zero
2571	transport_temp(:,:) = zero !tp_00 its a transport, not a temperature !!
2572	transport_temp_glo(:,:) = zero !tp_00
2573
2574	irrig_netereq(:) = zero
2575	irrig_needs(:,:) = zero
2576	irrig_actual(:,:) = zero
2577	irrig_deficit(:,:) = zero
2578	irrig_adduct(:,:) = zero
2579	totarea(:) = zero
2580	totflood(:) = zero
2581	!
2582	! Compute all the fluxes
2583	!
2584	DO ib=1,nbasmax
2585	DO ig=1,nbpt
2586	!
2587	totarea(ig) = totarea(ig) + routing_area(ig,ib)
2588	totflood(ig) = totflood(ig) + flood_reservoir(ig,ib)
2589	ENDDO
2590	ENDDO
2591	!
2592	!> The outflow fluxes from the three reservoirs are computed.
2593	!> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume.
2594	!> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid
2595	!> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula
2596	!> (Dingman, 1994; Ducharne et al., 2003).
2597	!> The resulting product of tcst (in s/km) and topo_resid (in km) represents the time constant (s)
2598	!> which is an e-folding time, the time necessary for the water amount
2599	!> in the stream reservoir to decrease by a factor e. Hence, it gives an order of
2600	!> magnitude of the travel time through this reservoir between
2601	!> the sub-basin considered and its downstream neighbor.
2602	!
2603	CALL groundwatertemp(nbpt, nbasmax, ngrnd, tempdiag, znt, dlt, fast_temp, slow_temp)
2604	!
2605	streamlimit(:) = zero
2606	!
2607	DO ib=1,nbasmax
2608	DO ig=1,nbpt
2609	IF ( route_tobasin(ig,ib) .GT. 0 ) THEN
2610	!
2611	! Each of the fluxes is limited by the water in the reservoir and a small margin
2612	! (min_reservoir) to avoid rounding errors.
2613	!
2614	flow = MIN(fast_reservoir(ig,ib)/(topo_resid(ig,ib)*fast_tcst/dt_routing),&
2615	& fast_reservoir(ig,ib)-min_sechiba)
2616	fast_flow(ig,ib) = MAX(flow, zero)
2617
2618	flow = MIN(slow_reservoir(ig,ib)/(topo_resid(ig,ib)*slow_tcst/dt_routing),&
2619	& slow_reservoir(ig,ib)-min_sechiba)
2620	slow_flow(ig,ib) = MAX(flow, zero)
2621
2622	! Need to adjust the reduction of the flow
2623	reduced = MAX(1-SQRT(MIN(flood_frac_bas(ig,ib),rfloodmax)), min_sechiba) ! Add the reduction flow parameter
2624	flow = stream_reservoir(ig,ib)/(stream_resid(ig,ib)stream_tcst/dt_routing)reduced
2625	flow = MIN(flow, stream_reservoir(ig,ib)-min_sechiba)
2626	stream_flow(ig,ib) = MAX(flow, zero)
2627	IF ( stream_flow(ig,ib) .GE. stream_reservoir(ig,ib)-min_sechiba .AND. stream_flow(ig,ib) > zero .AND. &
2628	& routing_area(ig,ib) > zero ) THEN
2629	streamlimit(ig) = streamlimit(ig)+1.0
2630	ENDIF
2631	!
2632	ELSE
2633	fast_flow(ig,ib) = zero
2634	slow_flow(ig,ib) = zero
2635	stream_flow(ig,ib) = zero
2636	ENDIF
2637	ENDDO
2638	ENDDO
2639	!-
2640	!- Compute the fluxes out of the floodplains and ponds if they exist.
2641	!-
2642	IF (do_floodplains .OR. doponds) THEN
2643	DO ig=1,nbpt
2644	IF (flood_frac(ig) .GT. min_sechiba) THEN
2645	!!!!
2646	! PONDS : not actualized
2647	!
2648	!flow = MIN(floodout(ig)totarea(ig)pond_frac(ig)/flood_frac(ig), pond_reservoir(ig)+totflood(ig))
2649	!pondex = MAX(flow - pond_reservoir(ig), zero)
2650	!pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex)
2651	!
2652	! If demand was over reservoir size, we will take it out from floodplains
2653	!
2654	!pond_excessflow(:) = zero
2655	!DO ib=1,nbasmax
2656	! pond_excessflow(ib) = MIN(pondex*flood_frac_bas(ig,ib)/(flood_frac(ig)-pond_frac(ig)),&
2657	! & flood_reservoir(ig,ib))
2658	! pondex = pondex - pond_excessflow(ib)
2659	!ENDDO
2660	!
2661	!IF ( pondex .GT. min_sechiba) THEN
2662	! WRITE(numout,*) "Unable to redistribute the excess pond outflow over the water available in the floodplain."
2663	! WRITE(numout,*) "Pondex = ", pondex
2664	! WRITE(numout,*) "pond_excessflow(:) = ", pond_excessflow(:)
2665	!ENDIF
2666	!
2667	DO ib=1,nbasmax
2668	!
2669	! when ponds actualized : add pond_excessflow to flow
2670	! This is the flow out of the reservoir due to ET (+ pond excessflow(ig), suppressed here)
2671	!flow = floodout(ig)routing_area(ig,ib)flood_frac_bas(ig,ib)/flood_frac(ig)
2672	flow = floodout(ig)routing_area(ig,ib)flood_frac_bas(ig,ib)
2673	!
2674	flood_reservoir(ig,ib) = flood_reservoir(ig,ib) - flow
2675	!
2676	!
2677	IF (flood_reservoir(ig,ib) .LT. min_sechiba) THEN
2678	flood_reservoir(ig,ib) = zero
2679	ENDIF
2680	IF (pond_reservoir(ig) .LT. min_sechiba) THEN
2681	pond_reservoir(ig) = zero
2682	ENDIF
2683	ENDDO
2684	ENDIF
2685	ENDDO
2686	ENDIF
2687
2688	!-
2689	!- Computing the drainage and outflow from floodplains
2690	!> Drainage from floodplains is depending on a averaged conductivity (k_litt)
2691	!> for saturated infiltration in the 'litter' layer. Flood_drainage will be
2692	!> a component of the total reinfiltration that leaves the routing scheme.
2693	!-
2694	IF (do_floodplains) THEN
2695	IF (dofloodinfilt) THEN
2696	DO ib=1,nbasmax
2697	DO ig=1,nbpt
2698	flood_drainage(ig,ib) = MAX(zero, MIN(flood_reservoir(ig,ib), &
2699	& flood_frac_bas(ig,ib)* routing_area(ig,ib) * k_litt(ig) * &
2700	& conduct_factor * dt_routing/one_day))
2701	flood_reservoir(ig,ib) = flood_reservoir(ig,ib) - flood_drainage(ig,ib)
2702	ENDDO
2703	ENDDO
2704	ELSE
2705	DO ib=1,nbasmax
2706	DO ig=1,nbpt
2707	flood_drainage(ig,ib) = zero
2708	ENDDO
2709	ENDDO
2710	ENDIF
2711	!> Outflow from floodplains is computed depending a delay. This delay is characterized by a time constant
2712	!> function of the surface of the floodplains and the product of topo_resid and flood_tcst. flood_tcst
2713	!> has been calibrated through observations in the Niger Inner Delta (D'Orgeval, 2006).
2714	!
2715	DO ib=1,nbasmax
2716	DO ig=1,nbpt
2717	IF ( route_tobasin(ig,ib) .GT. 0 ) THEN
2718	IF (flood_reservoir(ig,ib) .GT. min_sechiba) THEN
2719	flow = MIN(flood_reservoir(ig,ib)/(stream_resid(ig,ib)*flood_tcst/dt_routing),&
2720	& flood_reservoir(ig,ib)-min_sechiba)
2721	flow = MAX(flow, zero)
2722	ELSE
2723	flow = zero
2724	ENDIF
2725	flood_flow(ig,ib) = flow
2726	ELSE
2727	flood_flow(ig,ib) = zero
2728	ENDIF
2729	ENDDO
2730	ENDDO
2731	ELSE
2732	DO ib=1,nbasmax
2733	DO ig=1,nbpt
2734	flood_drainage(ig,ib) = zero
2735	flood_flow(ig,ib) = zero
2736	flood_reservoir(ig,ib) = zero
2737	ENDDO
2738	ENDDO
2739	ENDIF
2740
2741	!-
2742	!- Computing drainage and inflow for ponds
2743	!> Drainage from ponds is computed in the same way than for floodplains.
2744	!> Reinfiltrated fraction from the runoff (i.e. the outflow from the fast reservoir)
2745	!> is the inflow of the pond reservoir.
2746	!-
2747	IF (doponds) THEN
2748	! If used, the slope coef is not used in hydrol for water2infilt
2749	DO ib=1,nbasmax
2750	DO ig=1,nbpt
2751	pond_inflow(ig,ib) = fast_flow(ig,ib) * reinf_slope(ig)
2752	pond_drainage(ig,ib) = MIN(pond_reservoir(ig)*routing_area(ig,ib)/totarea(ig), &
2753	& pond_frac(ig)routing_area(ig,ib)k_litt(ig)*dt_routing/one_day)
2754	fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib)
2755	ENDDO
2756	ENDDO
2757	ELSE
2758	DO ib=1,nbasmax
2759	DO ig=1,nbpt
2760	pond_inflow(ig,ib) = zero
2761	pond_drainage(ig,ib) = zero
2762	pond_reservoir(ig) = zero
2763	ENDDO
2764	ENDDO
2765	ENDIF
2766
2767	source(:,:) = fast_flow(:,:) + slow_flow(:,:) + stream_flow(:,:)
2768	CALL downstreamsum(nbpt, nbasmax, source, transport)
2769	source(:,:) = fast_flow(:,:)fast_temp(:,:) + slow_flow(:,:)slow_temp(:,:) + &
2770	& stream_flow(:,:)*stream_temp(:,:)
2771	CALL downstreamsum(nbpt, nbasmax, source, transport_temp)
2772	!-
2773	!- Do the floodings - First initialize
2774	!-
2775	return_swamp(:,:)=zero
2776	floods(:,:)=zero
2777	!-
2778	!> Over swamp areas, a fraction of water (return_swamp) is withdrawn from the river depending on the
2779	!> parameter swamp_cst.
2780	!> It will be transferred into soil moisture and thus does not return directly to the river.
2781	!
2782	!- 1. Swamps: Take out water from the river to put it to the swamps
2783	!-
2784	!
2785	IF ( doswamps ) THEN
2786	tobeflooded(:) = swamp(:)
2787	DO ib=1,nbasmax
2788	DO ig=1,nbpt
2789	potflood(ig,ib) = transport(ig,ib)
2790	!
2791	IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > 0. .AND. floodtemp(ig) > tp_00 ) THEN
2792	!
2793	IF (routing_area(ig,ib) > tobeflooded(ig)) THEN
2794	floodindex = tobeflooded(ig) / routing_area(ig,ib)
2795	ELSE
2796	floodindex = 1.0
2797	ENDIF
2798	return_swamp(ig,ib) = swamp_cst * potflood(ig,ib) * floodindex
2799	!
2800	tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib)
2801	!
2802	ENDIF
2803	ENDDO
2804	ENDDO
2805	ENDIF
2806	!-
2807	!- 2. Floodplains: Update the reservoir with the flux computed above.
2808	!-
2809	IF ( do_floodplains ) THEN
2810	DO ig=1,nbpt
2811	DO ib=1,nbasmax
2812	IF (floodplains(ig, ib) .GT. min_sechiba .AND. floodtemp(ig) .GT. tp_00) THEN
2813	floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib)
2814	ENDIF
2815	ENDDO
2816	ENDDO
2817	ENDIF
2818	!
2819	! Update all reservoirs
2820	!> The slow and deep reservoir (slow_reservoir) collect the deep drainage whereas the
2821	!> fast_reservoir collects the computed surface runoff. Both discharge into a third reservoir
2822	!> (stream_reservoir) of the next sub-basin downstream.
2823	!> Water from the floodplains reservoir (flood_reservoir) flows also into the stream_reservoir of the next sub-basin downstream.
2824	!> Water that flows into the pond_reservoir is withdrawn from the fast_reservoir.
2825	!
2826	negslow = 0
2827	DO ig=1,nbpt
2828	DO ib=1,nbasmax
2829	!
2830	fast_reservoir(ig,ib) = fast_reservoir(ig,ib) + runoff(ig)*routing_area(ig,ib) - &
2831	& fast_flow(ig,ib) - pond_inflow(ig,ib)
2832	!
2833	slow_reservoir(ig,ib) = slow_reservoir(ig,ib) + drainage(ig)*routing_area(ig,ib) - &
2834	& slow_flow(ig,ib)
2835	!
2836	oldstream = stream_reservoir(ig, ib) * stream_temp(ig,ib)
2837	!
2838	stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_flow(ig,ib) + transport(ig,ib) - &
2839	& stream_flow(ig,ib) - return_swamp(ig,ib) - floods(ig,ib)
2840	!
2841	! Diagnostics of the stream reservoir
2842	!
2843	IF ( routing_area(ig,ib) > zero ) THEN
2844	! 1000 to transform kg into m^3
2845	htmp = stream_reservoir(ig,ib)*1000/routing_area(ig,ib)
2846	ewh(ig,ib) = 1.0/(1.0+htmp*hscale)
2847	ELSE
2848	ewh(ig,ib) = 1.0
2849	ENDIF
2850	!
2851	!reste du calcul
2852	!
2853	krelax = ewh(ig,ib)
2854	!
2855	den = 1.0/(1.0+dt_routing*krelax)
2856	IF ( stream_reservoir(ig,ib) > 1.e-6 ) THEN
2857	oldtemp = stream_temp(ig,ib)
2858	stream_temp(ig,ib) = den * dt_routing * krelax * fast_temp(ig,ib) + &
2859	& den * oldstream/stream_reservoir(ig,ib) + &
2860	& den * transport_temp(ig, ib)/stream_reservoir(ig,ib) - &
2861	& den * oldtemp*stream_flow(ig,ib)/stream_reservoir(ig,ib)
2862	!
2863	!Stream_temp [K], stream_reservoir [kg], WaterCp [J/g/K] yields tendencies in GJ/s
2864	!
2865	stemp_total_tend(ig,ib) = WaterCp1.e-6(stream_temp(ig,ib)*stream_reservoir(ig,ib) - oldstream)/dt_routing
2866	stemp_advec_tend(ig,ib) = WaterCp1.e-6(transport_temp(ig, ib) - oldtemp*stream_flow(ig,ib))/dt_routing
2867	stemp_relax_tend(ig,ib) = WaterCp1.e-6stream_reservoir(ig,ib)krelax(fast_temp(ig,ib)-stream_temp(ig,ib))
2868	ELSE
2869	stream_temp(ig,ib) = MAX(fast_temp(ig,ib), ZeroCelsius)
2870	stemp_total_tend(ig,ib) = zero
2871	stemp_advec_tend(ig,ib) = zero
2872	stemp_relax_tend(ig,ib) = zero
2873	ENDIF
2874	!
2875	flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + floods(ig,ib) - &
2876	& flood_flow(ig,ib)
2877	!
2878	pond_reservoir(ig) = pond_reservoir(ig) + pond_inflow(ig,ib) - pond_drainage(ig,ib)
2879	!
2880	IF ( flood_reservoir(ig,ib) .LT. zero ) THEN
2881	IF ( check_reservoir ) THEN
2882	WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ib, ". Problem is being corrected."
2883	WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir(ig,ib), floods(ig,ib), &
2884	& flood_flow(ig,ib)
2885	ENDIF
2886	stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_reservoir(ig,ib)
2887	flood_reservoir(ig,ib) = zero
2888	ENDIF
2889	!
2890	IF ( stream_reservoir(ig,ib) .LT. zero ) THEN
2891	IF ( check_reservoir ) THEN
2892	WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ib, ". Problem is being corrected."
2893	WRITE(numout,*) "stream_reservoir, flood_flow, transport : ", stream_reservoir(ig,ib), flood_flow(ig,ib), &
2894	& transport(ig,ib)
2895	WRITE(numout,*) "stream_flow, return_swamp, floods :", stream_flow(ig,ib), return_swamp(ig,ib), floods(ig,ib)
2896	ENDIF
2897	fast_reservoir(ig,ib) = fast_reservoir(ig,ib) + stream_reservoir(ig,ib)
2898	stream_reservoir(ig,ib) = zero
2899	ENDIF
2900	!
2901	IF ( fast_reservoir(ig,ib) .LT. zero ) THEN
2902	IF ( check_reservoir ) THEN
2903	WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ib, ". Problem is being corrected."
2904	WRITE(numout,*) "fast_reservoir, runoff, fast_flow, ponf_inflow : ", fast_reservoir(ig,ib), &
2905	&runoff(ig), fast_flow(ig,ib), pond_inflow(ig,ib)
2906	ENDIF
2907	slow_reservoir(ig,ib) = slow_reservoir(ig,ib) + fast_reservoir(ig,ib)
2908	fast_reservoir(ig,ib) = zero
2909	ENDIF
2910
2911	IF ( slow_reservoir(ig,ib) .LT. - min_sechiba ) THEN
2912	IF ( negslow < 20 ) THEN
2913	negslow = negslow + 1
2914	negig(negslow) = ig
2915	negib(negslow) = ib
2916	ENDIF
2917	ENDIF
2918
2919	ENDDO
2920	ENDDO
2921
2922	IF ( negslow > 0 ) THEN
2923	DO ier = 1,negslow
2924	ig = negig(ier)
2925	ib = negib(ier)
2926	WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig, ib,lalo(ig,:)
2927	WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', &
2928	& slow_reservoir(ig,ib), slow_flow(ig,ib), drainage(ig)
2929	WRITE(numout,*) 'WARNING : pondr, pond_inflow, pond_drainage', &
2930	& pond_reservoir(ig), pond_inflow(ig,ib), pond_drainage(ig,ib)
2931	CALL ipslerr_p(2, 'routing_hr_flow', 'WARNING negative slow_reservoir.','','')
2932	ENDDO
2933	ENDIF
2934
2935	totflood(:) = zero
2936	DO ig=1,nbpt
2937	DO ib=1,nbasmax
2938	totflood(ig) = totflood(ig) + flood_reservoir(ig,ib)
2939	ENDDO
2940	ENDDO
2941	!
2942	! ESTIMATE the flooded fraction
2943	!
2944	IF (do_floodplains .OR. doponds) THEN
2945	CALL routing_hr_flood(nbpt, flood_frac, totarea, totflood)
2946	ELSE
2947	flood_frac(:) = zero
2948	flood_height(:,:) = zero
2949	flood_frac_bas(:,:) = zero
2950	ENDIF
2951
2952
2953	!! ANTHONY : OVERFLOW
2954	!! CALCULATE TRANSFER BETWEEN FLOODPLAINS RESERVOIR
2955	IF (do_floodplains .AND. dofloodoverflow) Then
2956	! The overflow is repeated "overflow_repetition" times
2957	! This is in order to have more stability and
2958	! be able to use lower "overflow_tcst".
2959	DO ier = 1,overflow_repetition
2960	CALL routing_hr_overflow(nbpt, nbasmax)
2961	END DO
2962	! Once done we update the floodplains fraction and the floodplains height
2963	CALL routing_hr_flood(nbpt, flood_frac, totarea, totflood)
2964	END IF
2965
2966
2967	!-
2968	!- Compute the total reinfiltration and returnflow to the grid box
2969	!> A term of returnflow is computed including the water from the swamps that does not return directly to the river
2970	!> but will be put into soil moisture (see hydrol module).
2971	!> A term of reinfiltration is computed including the water that reinfiltrated from the ponds and floodplains areas.
2972	!> It will be put into soil moisture (see hydrol module).
2973	!-
2974	IF (do_floodplains .OR. doswamps .OR. doponds) THEN
2975	returnflow(:) = zero
2976	reinfiltration(:) = zero
2977	!
2978	DO ib=1,nbasmax
2979	DO ig=1,nbpt
2980	returnflow(ig) = returnflow(ig) + return_swamp(ig,ib)
2981	reinfiltration(ig) = reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib)
2982	ENDDO
2983	ENDDO
2984	!
2985	DO ig=1,nbpt
2986	returnflow(ig) = returnflow(ig)/totarea(ig)
2987	reinfiltration(ig) = reinfiltration(ig)/totarea(ig)
2988	ENDDO
2989	ELSE
2990	returnflow(:) = zero
2991	reinfiltration(:) = zero
2992	ENDIF
2993
2994	!
2995	! Compute the net irrigation requirement from Univ of Kassel
2996	!
2997	! This is a very low priority process and thus only applies if
2998	! there is some water left in the reservoirs after all other things.
2999	!
3000	!> The computation of the irrigation is performed here.
3001	!> * First step
3002	!> In a first time, the water requirements (irrig_netereq) by the crops for their optimal growth are calculated
3003	!> over each irrigated fraction (irrigated(ig)/totarea(ig)). It is the difference
3004	!> between the maximal water loss by the crops (transpot_mean) and the net water amount kept by the soil
3005	!> (precipitation and reinfiltration). Transpot_mean is computed in the routines enerbil and diffuco. It
3006	!> is derived from the effective transpiration parametrization under stress-free conditions, called potential transpiration.
3007	!> Crop_coef was used by a previous parametrization of irrigation in the code. Here, its value is equal to one.
3008	!> The crop coefficient was constant in space and time to represent a mean resistance of the vegetation to the potential evaporation.
3009	!> Now, the term crop_coef*Epot is substituted by transpot_mean (see Guimberteau et al., 2011).
3010	!> * Second step
3011	!> We compute irrigation needs in order to supply Irrig_netereq. Water for irrigation (irrig_actual) is withdrawn
3012	!> from the reservoirs. The amount of water is withdrawn in priority from the stream reservoir.
3013	!> If the irrigation requirement is higher than the water availability of the reservoir, water is withdrawn
3014	!> from the fast reservoir or, in the extreme case, from the slow reservoir.
3015	!> * Third step
3016	!> We compute a deficit in water for irrigation. If it is positive, irrigation (depending on water availibility in the reservoirs)
3017	!> has not supplied the crops requirements.
3018	!
3019	IF ( do_irrigation ) THEN
3020	DO ig=1,nbpt
3021	!
3022	IF ((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. &
3023	& (runoff(ig) .LT. min_sechiba) ) THEN
3024
3025	irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - &
3026	& (precip(ig)+reinfiltration(ig)) )
3027
3028	ENDIF
3029	!
3030	DO ib=1,nbasmax
3031	IF ( routing_area(ig,ib) .GT. 0 ) THEN
3032
3033	irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib)
3034
3035	irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),&
3036	& stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) )
3037
3038	slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + &
3039	& MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))))
3040
3041	fast_reservoir(ig,ib) = MAX( zero, &
3042	& fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))
3043
3044	stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) )
3045
3046	irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib)
3047
3048	ENDIF
3049	ENDDO
3050	!
3051	! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only).
3052	! If we find that then we create some adduction from that subbasin to the one where we need it for
3053	! irrigation.
3054	!
3055	!> If crops water requirements have not been supplied (irrig_deficit>0), we check if we cannot find the missing water
3056	!> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction
3057	!> from that subbasin to the one where we need it for irrigation.
3058	!>
3059	DO ib=1,nbasmax
3060
3061	stream_tot = SUM(stream_reservoir(ig,:))
3062
3063	DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)
3064
3065	fi = MAXLOC(stream_reservoir(ig,:))
3066	ib2 = fi(1)
3067
3068	irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), stream_reservoir(ig,ib2))
3069	stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)
3070	irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)
3071
3072	stream_tot = SUM(stream_reservoir(ig,:))
3073
3074	ENDDO
3075
3076	ENDDO
3077	!
3078	ENDDO
3079	!
3080	! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams
3081	! which can feed irrigation
3082	!
3083	!> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes
3084	!> to the one where we need it for irrigation.
3085	!
3086	IF (is_root_prc) THEN
3087	ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), &
3088	& irrig_adduct_glo(nbp_glo, nbasmax), stat=ier)
3089	ELSE
3090	ALLOCATE(irrig_deficit_glo(0, 0), stream_reservoir_glo(0, 0), &
3091	& irrig_adduct_glo(0, 0), stat=ier)
3092	ENDIF
3093	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_flow','Pb in allocate for irrig_deficit_glo, stream_reservoir_glo,...','','')
3094
3095	CALL gather(irrig_deficit, irrig_deficit_glo)
3096	CALL gather(stream_reservoir, stream_reservoir_glo)
3097	CALL gather(irrig_adduct, irrig_adduct_glo)
3098
3099	IF (is_root_prc) THEN
3100	!
3101	DO ig=1,nbp_glo
3102	! Only work if the grid box is smaller than 100x100km. Else the piplines we build
3103	! here would be too long to be reasonable.
3104	IF ( resolution_g(ig,1) < 100000. .AND. resolution_g(ig,2) < 100000. ) THEN
3105	DO ib=1,nbasmax
3106	!
3107	IF ( irrig_deficit_glo(ig,ib) > min_sechiba ) THEN
3108	!
3109	streams_around(:,:) = zero
3110	!
3111	DO in=1,NbNeighb
3112	ig2 = neighbours_g(ig,in)
3113	IF (ig2 .GT. 0 ) THEN
3114	streams_around(in,:) = stream_reservoir_glo(ig2,:)
3115	igrd(in) = ig2
3116	ENDIF
3117	ENDDO
3118	!
3119	IF ( MAXVAL(streams_around) .GT. zero ) THEN
3120	!
3121	ff=MAXLOC(streams_around)
3122	ig2=igrd(ff(1))
3123	ib2=ff(2)
3124	!
3125	IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. stream_reservoir_glo(ig2,ib2) > zero ) THEN
3126	adduction = MIN(irrig_deficit_glo(ig,ib), stream_reservoir_glo(ig2,ib2))
3127	stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction
3128	irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction
3129	irrig_adduct_glo(ig,ib) = irrig_adduct_glo(ig,ib) + adduction
3130	ENDIF
3131	!
3132	ENDIF
3133	!
3134	ENDIF
3135	!
3136	ENDDO
3137	ENDIF
3138	ENDDO
3139	!
3140	ENDIF
3141	!
3142
3143	CALL scatter(irrig_deficit_glo, irrig_deficit)
3144	CALL scatter(stream_reservoir_glo, stream_reservoir)
3145	CALL scatter(irrig_adduct_glo, irrig_adduct)
3146
3147	DEALLOCATE(irrig_deficit_glo, stream_reservoir_glo, irrig_adduct_glo)
3148
3149	ENDIF
3150
3151	!! Calculate the net water flow to each routing reservoir (in kg/dt)
3152	!! to further diagnose the corresponding water budget residu
3153	!! in routing_highres_main
3154
3155	netflow_fast_diag(:) = zero
3156	netflow_slow_diag(:) = zero
3157	netflow_stream_diag(:) = zero
3158
3159	DO ib=1,nbasmax
3160	DO ig=1,nbpt
3161	netflow_fast_diag(ig) = netflow_fast_diag(ig) + runoff(ig)*routing_area(ig,ib) &
3162	- fast_flow(ig,ib) - pond_inflow(ig,ib)
3163	netflow_slow_diag(ig) = netflow_slow_diag(ig) + drainage(ig)*routing_area(ig,ib) &
3164	- slow_flow(ig,ib)
3165	netflow_stream_diag(ig) = netflow_stream_diag(ig) + flood_flow(ig,ib) + transport(ig,ib) &
3166	- stream_flow(ig,ib) - return_swamp(ig,ib) - floods(ig,ib)
3167	ENDDO
3168	ENDDO
3169
3170	!! Grid cell averaging
3171	DO ig=1,nbpt
3172	netflow_fast_diag(ig) = netflow_fast_diag(ig)/totarea(ig)
3173	netflow_slow_diag(ig) = netflow_slow_diag(ig)/totarea(ig)
3174	netflow_stream_diag(ig) = netflow_stream_diag(ig)/totarea(ig)
3175	ENDDO
3176
3177	!
3178	!
3179	! Compute the fluxes which leave the routing scheme
3180	!
3181	! Lakeinflow is in Kg/dt
3182	! returnflow is in Kg/m^2/dt
3183	!
3184	hydrographs(:) = zero
3185	hydrotemp(:) = zero
3186	HTUhgmon(:,:) = zero
3187	HTUtempmon(:,:) = zero
3188	slowflow_diag(:) = zero
3189	fast_diag(:) = zero
3190	slow_diag(:) = zero
3191	stream_diag(:) = zero
3192	flood_diag(:) = zero
3193	pond_diag(:) = zero
3194	irrigation(:) = zero
3195	!
3196	!
3197	DO ib=1,nbasmax
3198	!
3199	DO ig=1,nbpt
3200	!
3201	DO im=1,nbasmon
3202	IF (HTUdiag_loc(ig,im) > 0 .AND. HTUdiag_loc(ig,im) .EQ. ib ) THEN
3203	HTUhgmon(ig,im) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib)
3204	HTUtempmon(ig,im) = stream_temp(ig,ib)
3205	ENDIF
3206	ENDDO
3207	!
3208	IF (hydrodiag(ig) == ib) THEN
3209	hydrographs(ig) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib)
3210	hydrotemp(ig) = stream_temp(ig,ib)
3211	slowflow_diag(ig) = slowflow_diag(ig) + slow_flow(ig,ib)
3212	ENDIF
3213	fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib)
3214	slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib)
3215	stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib)
3216	flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib)
3217	irrigation (ig) = irrigation (ig) + irrig_actual(ig,ib) + irrig_adduct(ig,ib)
3218	ENDDO
3219	ENDDO
3220	!
3221	DO ig=1,nbpt
3222	fast_diag(ig) = fast_diag(ig)/totarea(ig)
3223	slow_diag(ig) = slow_diag(ig)/totarea(ig)
3224	stream_diag(ig) = stream_diag(ig)/totarea(ig)
3225	flood_diag(ig) = flood_diag(ig)/totarea(ig)
3226	pond_diag(ig) = pond_reservoir(ig)/totarea(ig)
3227	!
3228	irrigation(ig) = irrigation(ig)/totarea(ig)
3229	!
3230	! The three output types for the routing : endoheric basins,, rivers and
3231	! diffuse coastal flow.
3232	!
3233	lakeinflow(ig) = transport(ig,nbasmax+1)
3234	coastalflow(ig) = transport(ig,nbasmax+2)
3235	riverflow(ig) = transport(ig,nbasmax+3)
3236	!
3237	ENDDO
3238	!
3239	flood_res = flood_diag + pond_diag
3240
3241
3242	!! Remove water from lake reservoir if it exceeds the maximum limit and distribute it
3243	!! uniformly over all possible the coastflow gridcells
3244
3245	! Calculate lake_overflow and remove it from lake_reservoir
3246	DO ig=1,nbpt
3247	lake_overflow(ig) = MAX(0., lake_reservoir(ig) - max_lake_reservoir*totarea(ig))
3248	lake_reservoir(ig) = lake_reservoir(ig) - lake_overflow(ig)
3249	END DO
3250	! Transform lake_overflow from kg/grid-cell/dt_routing into kg/m^2/s
3251	CALL xios_orchidee_send_field("lake_overflow",lake_overflow(:)/totarea(:)/dt_routing)
3252
3253	! Calculate the sum of the lake_overflow and distribute it uniformly over all gridboxes
3254	CALL gather(lake_overflow,lake_overflow_g)
3255	IF (is_root_prc) THEN
3256	total_lake_overflow=SUM(lake_overflow_g)
3257	END IF
3258	CALL bcast(total_lake_overflow)
3259
3260	! Distribute the lake_overflow uniformly over all coastal gridcells
3261	! lake_overflow_coast is only calculated to be used as diagnostics if needed
3262	DO ig=1,nbpt
3263	coastalflow(ig) = coastalflow(ig) + total_lake_overflow/nb_coast_gridcells * mask_coast(ig)
3264	lake_overflow_coast(ig) = total_lake_overflow/nb_coast_gridcells * mask_coast(ig)
3265	END DO
3266	! Transform from kg/grid-cell/dt_routing into m^3/grid-cell/s to match output unit of coastalflow
3267	CALL xios_orchidee_send_field("lake_overflow_coast",lake_overflow_coast/mille/dt_routing)
3268
3269
3270	END SUBROUTINE routing_hr_flow
3271	!
3272	!! ================================================================================================================================
3273	!! SUBROUTINE : groundwatertemp
3274	!!
3275	!>\BRIEF : This subroutine computes the temperature of the groundwater leaving the HTU
3276	!!
3277	!! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir
3278	!! is based on a maximum lake evaporation rate (maxevap_lake). \n
3279	!!
3280	!! RECENT CHANGE(S): None
3281	!!
3282	!! MAIN OUTPUT VARIABLE(S):
3283	!!
3284	!! REFERENCES : None
3285	!!
3286	!! FLOWCHART :None
3287	!! \n
3288	!_ ================================================================================================================================
3289	!-
3290	SUBROUTINE groundwatertemp(nbpt, nbasmax, nl, tempdiag, lev, dlz, fast_temp, slow_temp)
3291	! INPUT
3292	INTEGER(i_std), INTENT(in) :: nbpt, nbasmax, nl
3293	REAL(r_std), INTENT(in) :: tempdiag(nbpt,nl)
3294	REAL(r_std), INTENT(in) :: lev(nl), dlz(nl)
3295	REAL(r_std), INTENT(inout) :: slow_temp(nbpt,nbasmax), fast_temp(nbpt,nbasmax)
3296	! OUTPUT
3297	! LOCAL
3298	INTEGER(i_std) :: ig, ib, im
3299	REAL(r_std) :: sw
3300	REAL(r_std) :: rw(nl), dw(nl)
3301	LOGICAL, SAVE :: alltop=.FALSE.
3302	LOGICAL, SAVE :: FirstCall=.TRUE.
3303	!
3304	IF ( FirstCall ) THEN
3305	!Config Key = ROUTING_ALLTOPT
3306	!Config Desc = Should drainage have the temperature of the top soil (0.3m) ?
3307	!Config Def = False
3308	!Config Help = The default behaviour of the scheme is that runoff has the temperature
3309	!Config Help of the top 30 cm of soil. Drainage will have the temperature of the lowest
3310	!Config Help soil layer (3-17m). If set to True this flag will give drainage the same
3311	!Config Help temperature as runoff.
3312	!Config Units = Logical
3313	alltop=.FALSE.
3314	CALL getin_p('ROUTING_ALLTOPT', alltop)
3315	!
3316	WRITE(numout,*) "Runoff will have the average soil temperature of layers from ", runofftempdepth(1),&
3317	& " to ", runofftempdepth(2), "[m]"
3318	!
3319	IF ( alltop ) THEN
3320	WRITE(numout,*) "Drainage will have the average soil temperature of layers from ", runofftempdepth(1),&
3321	& " to ", runofftempdepth(2), "[m]"
3322	ELSE
3323	WRITE(numout,*) "Drainage will have the average soil temperature of layers from ", drainagetempdepth(1),&
3324	& " to ", MIN(drainagetempdepth(2), SUM(dlz)), "[m]"
3325	ENDIF
3326	FirstCall=.FALSE.
3327	ENDIF
3328	!
3329	CALL tempdepthweight(nl, dlz, runofftempdepth(1), runofftempdepth(2), rw)
3330	CALL tempdepthweight(nl, dlz, drainagetempdepth(1), MIN(drainagetempdepth(2), SUM(dlz)), dw)
3331	!
3332	slow_temp(:,:) = zero
3333	fast_temp(:,:) = zero
3334	! Compute for each HTU the temperature of runoff and drainage water.
3335	DO im = 1,nl
3336	DO ib=1,nbasmax
3337	DO ig=1,nbpt
3338	fast_temp(ig,ib) = fast_temp(ig,ib) + tempdiag(ig,im)*rw(im)
3339	! The option to have drainage water at the same temperature as runoff
3340	IF ( alltop ) THEN
3341	slow_temp(ig,ib) = slow_temp(ig,ib) + tempdiag(ig,im)*rw(im)
3342	ELSE
3343	slow_temp(ig,ib) = slow_temp(ig,ib) + tempdiag(ig,im)*dw(im)
3344	ENDIF
3345	ENDDO
3346	ENDDO
3347	ENDDO
3348
3349	END SUBROUTINE groundwatertemp
3350
3351	SUBROUTINE tempdepthweight(n, dz, top, bot, w)
3352	! Input
3353	INTEGER(i_std), INTENT(in) :: n
3354	REAL(r_std), INTENT(in) :: dz(n)
3355	REAL(r_std), INTENT(in) :: top, bot
3356	! Output
3357	REAL(r_std), INTENT(out) :: w(n)
3358	! Local
3359	INTEGER(i_std) :: i
3360	REAL(r_std) :: sw
3361	w(:) = zero
3362	sw = zero
3363	DO i=1,n
3364	w(i) = MAX(zero, MIN(sw+dz(i), bot) - MAX(top, sw))
3365	sw = sw + dz(i)
3366	ENDDO
3367	w(:) = w(:)/(bot-top)
3368	END SUBROUTINE tempdepthweight
3369
3370	!! ================================================================================================================================
3371	!! SUBROUTINE : downstreamsum
3372	!!
3373	!>\BRIEF : This subroutine sums the input variables onto the downstream HTU in the river graph.
3374	!!
3375	!! DESCRIPTION : We assume that the downstream HTU is defined by route_togrid and route_tobas. As these
3376	!! donwstream HTU can be on another processor we do this job on the root processor. So before we need to
3377	!! transfer all the data onto that processor and then redistribute the result.
3378	!! Keep in mind that if an HTU does not exit then route_tobas = 0. So the result array needs
3379	!! to have this index. The end of the rivers are between nbmax+1 and nbmax+3 so this indexing space is also
3380	!! needed in the result array.
3381	!!
3382	!! RECENT CHANGE(S): None
3383	!!
3384	!! MAIN OUTPUT VARIABLE(S):
3385	!!
3386	!! REFERENCES : None
3387	!!
3388	!! FLOWCHART :None
3389	!! \n
3390	!_ ================================================================================================================================
3391	!-
3392	SUBROUTINE downstreamsum(nbpt, nbmax, v, t)
3393	! Input
3394	INTEGER(i_std), INTENT(in) :: nbpt, nbmax
3395	REAL(r_std), INTENT(in), DIMENSION(nbpt, nbmax) :: v
3396	! Output
3397	REAL(r_std), INTENT(out), DIMENSION(nbpt, 0:nbmax+3) :: t
3398	!
3399	! Local
3400	!
3401	INTEGER(i_std) :: ig, ib, rtg, rtb
3402	INTEGER(i_std) :: ier
3403	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: v_g, t_g
3404	!
3405	! Allocate memory if needed. Should only happen only once in order to reduce computing time.
3406	!
3407	IF ( .NOT. ALLOCATED(v_g) ) THEN
3408	IF (is_root_prc) THEN
3409	ALLOCATE(v_g(nbp_glo,nbmax), stat=ier)
3410	IF (ier /= 0) CALL ipslerr_p(3,'downstreamsum','Pb in allocate for v_g','','')
3411	ELSE
3412	ALLOCATE(v_g(1,1))
3413	ENDIF
3414	ENDIF
3415	IF ( .NOT. ALLOCATED(t_g) ) THEN
3416	IF (is_root_prc) THEN
3417	ALLOCATE(t_g(nbp_glo,0:nbmax+3), stat=ier)
3418	IF (ier /= 0) CALL ipslerr_p(3,'downstreamsum','Pb in allocate for t_g','','')
3419	ELSE
3420	ALLOCATE(t_g(1,1))
3421	ENDIF
3422	ENDIF
3423	!
3424	! Gather the source variable on the root processor.
3425	!
3426	CALL gather(v, v_g)
3427	!
3428	! The downstream sum is performed only on the root processor.
3429	!
3430	IF (is_root_prc) THEN
3431	t_g(:,:) = zero
3432	DO ib=1,nbmax
3433	DO ig=1,nbp_glo
3434	rtg = route_togrid_glo(ig,ib)
3435	rtb = route_tobasin_glo(ig,ib)
3436	t_g(rtg,rtb) = t_g(rtg,rtb) + v_g(ig,ib)
3437	ENDDO
3438	ENDDO
3439	ENDIF
3440	!
3441	! Redistribute the downstream field to the all processors.
3442	!
3443	CALL scatter(t_g, t)
3444	!
3445	END SUBROUTINE downstreamsum
3446	!! ================================================================================================================================
3447	!! SUBROUTINE : routing_hr_flood
3448	!!
3449	!>\BRIEF : This subroutine estimate the flood fraction and the flood height for each HTU
3450	!!
3451	!! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir
3452	!! is based on a maximum lake evaporation rate (maxevap_lake). \n
3453	!!
3454	!! RECENT CHANGE(S): None
3455	!!
3456	!! MAIN OUTPUT VARIABLE(S):
3457	!!
3458	!! REFERENCES : None
3459	!!
3460	!! FLOWCHART :None
3461	!! \n
3462	!_ ================================================================================================================================
3463	!-
3464	SUBROUTINE routing_hr_flood(nbpt, flood_frac, totarea, totflood)
3465	!
3466	IMPLICIT NONE
3467	!
3468	!! INPUT VARIABLES
3469	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
3470	REAL(r_std), INTENT(in), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg)
3471	REAL(r_std), INTENT(in), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2)
3472	!! Flooded fraction of the grid box (unitless;0-1)
3473	!
3474	!! OUTPUT VARIABLES
3475	REAL(r_std), INTENT(inout) :: flood_frac(nbpt)
3476
3477	!
3478	!! LOCAL VARIABLES
3479	INTEGER(i_std) :: ig, ib !! Indices (unitless)
3480	REAL(r_std) :: diff, voltemp !! Discharge reduction due to floodplains
3481	!_ ================================================================================================================================
3482	!
3483	!
3484	! Initialize the variables
3485	flood_frac(:) = zero
3486	flood_height(:,:) = zero
3487	flood_frac_bas(:,:) = zero
3488	DO ig=1, nbpt
3489	IF (totflood(ig) .GT. min_sechiba) THEN
3490	DO ib=1,nbasmax
3491	IF (floodplains(ig,ib) .GT. min_sechiba) THEN
3492	! We have to convert h0 to m and the flood_reservoir in m^3
3493	flood_frac_bas(ig,ib) = ((fp_beta(ig,ib)+un) * flood_reservoir(ig,ib) / 1000) / ( floodcri(ig,ib) / 1000 * floodplains(ig,ib))
3494	flood_frac_bas(ig,ib) = (flood_frac_bas(ig,ib)) ** (fp_beta(ig,ib)/(fp_beta(ig,ib)+1))
3495	flood_frac_bas(ig,ib) = MIN(flood_frac_bas(ig,ib), floodplains(ig,ib)/ routing_area(ig,ib) )
3496
3497	! flood_height is in mm
3498	! there is two cases: flood_height < h0, flood_height >= h0 (this corresponds to flood_frac_bas = 1 )
3499	IF ( flood_frac_bas(ig,ib) .EQ. floodplains(ig,ib) / routing_area(ig,ib) ) THEN
3500	! voltemp is on m^3
3501	! Calculation of volume corresponding to h0
3502	voltemp = floodplains(ig,ib)/(fp_beta(ig,ib)+un) * ( floodcri(ig,ib) / 1000 )
3503	voltemp = flood_reservoir(ig,ib) / 1000 - voltemp
3504	! flood height is in mm
3505	flood_height(ig, ib) = voltemp / floodplains(ig,ib) * 1000 + floodcri(ig,ib)
3506	ELSE
3507	! flood height is in mm
3508	flood_height(ig, ib) = (flood_frac_bas(ig,ib)) ** (1/fp_beta(ig,ib)) * floodcri(ig,ib)
3509	END IF
3510	ENDIF
3511	ENDDO
3512	ENDIF
3513
3514	DO ib=1,nbasmax
3515	flood_frac(ig) = flood_frac(ig) + flood_frac_bas(ig,ib) * routing_area(ig,ib) / totarea(ig)
3516	END DO
3517	flood_frac(ig) = flood_frac(ig) + pond_frac(ig)
3518	!
3519	ENDDO
3520
3521	END SUBROUTINE routing_hr_flood
3522	!
3523	!! ================================================================================================================================
3524	!! SUBROUTINE : routing_hr_overflow
3525	!!
3526	!>\BRIEF : This subroutine performs the overflow fluxes
3527	!!
3528	!! DESCRIPTION (definitions, functional, design, flags): \n
3529	!!
3530	!! RECENT CHANGE(S): None
3531	!!
3532	!! MAIN OUTPUT VARIABLE(S):
3533	!!
3534	!! REFERENCES : None
3535	!!
3536	!! FLOWCHART :None
3537	!! \n
3538	!_ ================================================================================================================================
3539	!-
3540	SUBROUTINE routing_hr_overflow(nbpt, nbasmax)
3541	!
3542	IMPLICIT NONE
3543	!
3544	!! INPUT VARIABLES
3545	INTEGER(i_std), INTENT(in) :: nbpt,nbasmax !! Domain size (unitless)
3546	!
3547	!! LOCAL VARIABLES
3548	REAL(r_std), DIMENSION(nbpt,nbasmax) :: transport_overflow !! Water transport between floodplains - flood overflow (kg/dt)
3549	REAL(r_std), DIMENSION(nbp_glo,nbasmax) :: transport_overflow_glo !! Water transport between floodplains - flood overflow (kg/dt)
3550	REAL(r_std), DIMENSION(nbpt,nbasmax) :: overflow_loss !! Water loss from flood overflow (kg/dt)
3551	REAL(r_std), DIMENSION(nbp_glo,nbasmax) :: overflow_loss_glo !! Water loss from flood overflow (kg/dt)
3552	!
3553	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_height_g !! Floodplains height (m)
3554	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas_g !! Fraction of the HTU flooded
3555	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir_glo !! Water amount in the stream reservoir (kg)
3556	!
3557	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: DH,DH_temp !! Difference of height - flood overflow (kg/dt)
3558	!
3559	INTEGER(i_std) :: numflood !!
3560	!
3561	INTEGER(i_std) :: ig, ib, inf,inb,ing !! Indices (unitless)
3562	REAL(r_std) :: diff !! Discharge reduction due to floodplains
3563	REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt)
3564	REAL(r_std) :: dorog !! Discharge reduction due to floodplains
3565	INTEGER(i_std) :: ier !! Error handling
3566
3567
3568	!_ ================================================================================================================================
3569	!
3570	!! ANTHONY : OVERFLOW
3571	!! CALCULATE TRANSFER BETWEEN FLOODPLAINS RESERVOIR
3572	IF (is_root_prc) THEN
3573	ALLOCATE( flood_height_g(nbp_glo, nbasmax), flood_frac_bas_g(nbp_glo, nbasmax), stat=ier)
3574	ALLOCATE( flood_reservoir_glo(nbp_glo, nbasmax), stat=ier)
3575	ELSE
3576	ALLOCATE( flood_height_g(1,1), flood_frac_bas_g(1,1), stat=ier)
3577	ALLOCATE( flood_reservoir_glo(1, 1), stat=ier)
3578	ENDIF
3579	!
3580	IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_flow','Pb in allocate for flood_height_glo/floog_frac_glo','','')
3581	!
3582	CALL gather(flood_height,flood_height_g)
3583	CALL gather(flood_frac_bas,flood_frac_bas_g)
3584	CALL gather(flood_reservoir,flood_reservoir_glo)
3585	!
3586	IF (is_root_prc) THEN
3587	transport_overflow_glo(:,:) = 0
3588	overflow_loss_glo(:,:) = 0
3589	DO ib=1,nbasmax
3590	DO ig=1,nbp_glo
3591	IF ( floodplains_glo(ig,ib)/routing_area_glo(ig,ib) .GT. 0.5) THEN
3592	numflood = 0 ! Number of inflows for overflow
3593	ALLOCATE(DH(route_innum_glo(ig,ib)))
3594	DH(:) = 0
3595	DH_temp(:) = -1
3596	DO inf=1,route_innum_glo(ig,ib)
3597	ing = route_ingrid_glo(ig,ib,inf)
3598	inb = route_inbasin_glo(ig,ib,inf)
3599	IF ( floodplains_glo(ing,inb)/routing_area_glo(ing,inb) .GT. 0 ) THEN
3600	! Minimum of deltaorog is defined at lim_floodcri (0.3 m
3601	! can be used).
3602	dorog = MAX(orog_min_glo(ing,inb)- orog_min_glo(ig,ib), lim_floodcri)
3603	! flood_height is in mm and orog min in m
3604	diff = (flood_height_g(ig,ib)- flood_height_g(ing,inb))/1000 - dorog
3605	DH(inf) = max(diff, 0.)
3606	!
3607	! Flux is estimated via floodplains_glo
3608	! Then factor 1000 is to convert m^3 to kg
3609	! OVERFLOW_TCST is in seconds
3610	flow = DH(inf) * (floodplains_glo(ig,ib)* floodplains_glo(ing,inb))/(floodplains_glo(ig,ib)+floodplains_glo(ing,inb))1000 / overflow_tcst dt_routing / one_day
3611	transport_overflow_glo(ing,inb) = transport_overflow_glo(ing,inb) + flow
3612	overflow_loss_glo(ig,ib) = overflow_loss_glo(ig,ib) + flow
3613	END IF
3614	END DO
3615	DEALLOCATE(DH)
3616	END IF
3617	ENDDO
3618	ENDDO
3619	END IF
3620	! Send to local variables
3621	CALL scatter(transport_overflow_glo, transport_overflow)
3622	CALL scatter(overflow_loss_glo, overflow_loss)
3623	! Apply the volume changes
3624	DO ig=1,nbpt
3625	DO ib=1,nbasmax
3626	IF ( floodplains(ig,ib) .GT. 0 ) THEN
3627	flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + transport_overflow(ig,ib) - overflow_loss(ig,ib)
3628	! NEED to check if flood reservoir is less than 0, this may be a critical issue
3629	! Solved by an adequate use of an higher overflow time constant
3630	! To obtain the same result as with a lower overflow parameter
3631	! -> repeat a few time the operation with and higher overflow parameter
3632	IF ( flood_reservoir(ig,ib) .LT. 0 ) THEN
3633
3634	WRITE(,) "Issue of flood reservoir < 0 due to overflow at ", ig, ib
3635	stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + stream_reservoir(ig,ib) ! + because negative !
3636	flood_reservoir(ig,ib) = 0
3637	END IF
3638	END IF
3639	END DO
3640	END DO
3641	DEALLOCATE( flood_height_g, flood_frac_bas_g)
3642
3643	END SUBROUTINE routing_hr_overflow
3644	!
3645	!! ================================================================================================================================
3646	!! SUBROUTINE : routing_hr_lake
3647	!!
3648	!>\BRIEF : This subroutine stores water in lakes so that it does not cycle through the runoff.
3649	!! For the moment it only works for endoheric lakes but I can be extended in the future.
3650	!!
3651	!! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir
3652	!! is based on a maximum lake evaporation rate (maxevap_lake). \n
3653	!!
3654	!! RECENT CHANGE(S): None
3655	!!
3656	!! MAIN OUTPUT VARIABLE(S):
3657	!!
3658	!! REFERENCES : None
3659	!!
3660	!! FLOWCHART :None
3661	!! \n
3662	!_ ================================================================================================================================
3663
3664	SUBROUTINE routing_hr_lake(nbpt, dt_routing, lakeinflow, humrel, return_lakes)
3665	!
3666	IMPLICIT NONE
3667	!
3668	!! INPUT VARIABLES
3669	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
3670	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
3671	REAL(r_std), INTENT(out) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt)
3672	REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless)
3673	!
3674	!! OUTPUT VARIABLES
3675	REAL(r_std), INTENT(out) :: return_lakes(nbpt) !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
3676	!
3677	!! LOCAL VARIABLES
3678	INTEGER(i_std) :: ig !! Indices (unitless)
3679	REAL(r_std) :: refill !!
3680	REAL(r_std) :: total_area !! Sum of all the surfaces of the basins (m^2)
3681
3682	!_ ================================================================================================================================
3683	!
3684	!
3685	DO ig=1,nbpt
3686	!
3687	total_area = SUM(routing_area(ig,:))
3688	!
3689	lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig)
3690
3691	IF ( doswamps ) THEN
3692	! Calculate a return flow that will be extracted from the lake reservoir and reinserted in the soil in hydrol
3693	! Uptake in Kg/dt
3694	refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area)
3695	return_lakes(ig) = MIN(refill, lake_reservoir(ig))
3696	lake_reservoir(ig) = lake_reservoir(ig) - return_lakes(ig)
3697	! Return in Kg/m^2/dt
3698	return_lakes(ig) = return_lakes(ig)/total_area
3699	ELSE
3700	return_lakes(ig) = zero
3701	ENDIF
3702
3703	! This is the volume of the lake scaled to the entire grid.
3704	! It would be better to scale it to the size of the lake
3705	! but this information is not yet available.
3706	lake_diag(ig) = lake_reservoir(ig)/total_area
3707
3708	lakeinflow(ig) = lakeinflow(ig)/total_area
3709
3710	ENDDO
3711	!
3712	END SUBROUTINE routing_hr_lake
3713	!
3714	!! ================================================================================================================================
3715	!! SUBROUTINE : routing_hr_basins_p
3716	!!
3717	!>\BRIEF This routing read the file created by RoutingPreProc : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp
3718	!!
3719	!! DESCRIPTION (definitions, functional, design, flags) : None
3720	!! Once the atmospheric grid is defined and the land/sea mask set, RoutingPreProc has to used to generate the
3721	!! HTU graphs for the domain. This can be done either on the basis of the HydroSHEDS, MERIT or the old VÃ¶rÃ¶smarty map
3722	!! of catchments. During this step all the information will be created to allow ORCHIDEE to route the water and
3723	!! and monitor the flows at given stations.
3724	!! For the moment the ROUTING_FILE (Perhaps to renamed RoutingGraph) is read using IOIPSL but that should evolve toward XIOS.
3725	!!
3726	!! RECENT CHANGE(S): None
3727	!!
3728	!! MAIN OUTPUT VARIABLE(S):
3729	!!
3730	!! REFERENCES : None
3731	!!
3732	!! FLOWCHART : None
3733	!! \n
3734	!_ ================================================================================================================================
3735
3736	SUBROUTINE routing_hr_basins_p(nbpt, lalo, neighbours, resolution, contfrac)
3737	!
3738	IMPLICIT NONE
3739	!
3740	!! INPUT VARIABLES
3741	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
3742	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
3743	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb) !! Vector of neighbours for each grid point (1=North and then clockwise) (unitless)
3744	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m)
3745	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
3746	!
3747	! LOCAL
3748	!
3749	INTEGER(i_std) :: iml, jml, lml, tml
3750	INTEGER(i_std) :: i, j, ni, fid, ib, ig, ic, ign, ibn, og, ob, ier, im
3751	REAL(r_std) :: corr
3752	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: tmpvar_glo
3753	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmpvar
3754	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon, lat, landindex
3755	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: indextab
3756	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: landfileindex
3757	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: land2land
3758	INTEGER(i_std) :: nbhtumon
3759	!
3760	!_ ================================================================================================================================
3761	!
3762	!
3763	!
3764	IF (is_root_prc) THEN
3765	!
3766	CALL flininfo(graphfilename, iml, jml, lml, tml, fid)
3767	!
3768	IF (iml .NE. iim_g .AND. jml .NE. jjm_g ) THEN
3769	CALL ipslerr(3,'routing_hr_basins_p',&
3770	'The routing graph file does not have the right dimensions for the model.', &
3771	'Are you sure you are using the right routing graph file ?', ' ')
3772	ENDIF
3773	!
3774	!
3775	ALLOCATE(tmpvar_glo(iml,jml,nbasmax))
3776	ALLOCATE(tmpvar(iml,jml))
3777	ALLOCATE(lon(iml,jml))
3778	ALLOCATE(lat(iml,jml))
3779	ALLOCATE(landindex(iml,jml))
3780	ALLOCATE(indextab(iml,jml))
3781	ALLOCATE(landfileindex(iml,jml))
3782	ALLOCATE(land2land(iml*jml))
3783	!
3784	CALL flinget(fid, 'lon', iml, jml, 1, tml, 1, 0, lon)
3785	CALL flinget(fid, 'lat', iml, jml, 1, tml, 1, 0, lat)
3786	CALL flinget(fid, 'nbpt_glo', iml, jml, 1, tml, 1, 0, landindex)
3787	!
3788	! Replace NaN and other undef values
3789	!
3790	DO i=1,iml
3791	DO j=1,jml
3792	IF ( landindex(i,j) /= landindex(i,j) .OR. landindex(i,j) >= undef_graphfile) THEN
3793	landindex(i,j) = -1
3794	ENDIF
3795	ENDDO
3796	ENDDO
3797	!
3798	! Compute land index for file data. Information could be in file !
3799	!
3800	ni=NINT(MAXVAL(landindex))
3801	IF ( ni .NE. nbp_glo) THEN
3802	WRITE(numout,*) "Error routing_hr_basins_p : ni, nbp_glo : ", ni, nbp_glo, undef_graphfile
3803	CALL ipslerr(3,'routing_hr_basins_p',&
3804	'The routing graph file does not have the same number', &
3805	'of land points as the model.',&
3806	' ')
3807	ENDIF
3808	!
3809	CALL routing_hr_indexfilegrid(iml, jml, nbp_glo, lon, lat, landindex, indextab, land2land)
3810	!
3811	CALL flinget(fid, 'basin_area', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3812	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, routing_area_glo, &
3813	& zero)
3814
3815	IF ( do_floodplains ) THEN
3816	CALL flinget(fid, 'basin_floodp', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3817	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, floodplains_glo, &
3818	& zero)
3819	!
3820	CALL flinget(fid, 'floodcri', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3821	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, floodcri_glo, &
3822	& un)
3823	!
3824	CALL flinget(fid, 'basin_beta_fp', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3825	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, fp_beta_glo, &
3826	& un)
3827	END IF
3828
3829	CALL flinget(fid, 'topoindex', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3830	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, topo_resid_glo, &
3831	& undef_graphfile)
3832
3833	IF ( graphfile_version >= 2.0) THEN
3834	CALL flinget(fid, 'topoindex_stream', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3835	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, stream_resid_glo, &
3836	& undef_graphfile)
3837	CALL ipslerr(1,'routing_hr_basins_p',&
3838	'The topoindex_stream variable was found in routing_graph.nc', &
3839	'It will be used the topographic index of the stream store.',&
3840	' ')
3841	ELSE
3842	stream_resid_glo(:,:) = topo_resid_glo(:,:)
3843	ENDIF
3844	stream_maxresid=MAXVAL(stream_resid_glo, MASK=stream_resid_glo .LT. undef_graphfile)
3845
3846	CALL flinget(fid, 'basinid', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3847	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, global_basinid_glo, &
3848	& undef_int)
3849
3850	CALL flinget(fid, 'routetogrid', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3851	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, route_togrid_glo, &
3852	& undef_int)
3853	CALL routing_hr_convertlandpts(nbp_glo, nbasmax, land2land, route_togrid_glo)
3854
3855	CALL flinget(fid, 'routetobasin', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3856	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, route_tobasin_glo, 0)
3857
3858	CALL flinget(fid, 'routenbintobas', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3859	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, route_nbintobas_glo, 0)
3860
3861	!!
3862	IF ( dofloodoverflow ) THEN
3863	CALL flinget(fid, 'basin_orog_min', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo)
3864	CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, orog_min_glo, un)
3865	END IF
3866	!!
3867	IF ( graphfile_version >= 2.6) THEN
3868	CALL flinget(fid, 'gridrephtu', iml, jml, 1, tml, 1, 0, tmpvar)
3869	CALL routing_hr_landgather(iml, jml, nbp_glo, indextab, tmpvar, hydrodiag_glo, -1)
3870	ELSE
3871	hydrodiag_glo(:) = 1
3872	ENDIF
3873	!!
3874	IF ( MonitoringinGraph ) THEN
3875	CALL flinget(fid, 'HTUmonitor', iml, jml, nbasmon, tml, 1, 0, tmpvar_glo)
3876	CALL routing_hr_landgather(iml, jml, nbasmon, nbp_glo, indextab, tmpvar_glo, HTUdiag_glo, -1)
3877	ELSE
3878	HTUdiag_glo(:,:) = -1
3879	ENDIF
3880	!
3881	CALL flinclo(fid)
3882	DEALLOCATE(indextab)
3883	DEALLOCATE(lon)
3884	DEALLOCATE(lat)
3885	DEALLOCATE(tmpvar_glo)
3886	DEALLOCATE(tmpvar)
3887	!
3888	! Convert floodplains fraction into floodplains surface
3889	IF ( do_floodplains ) THEN
3890	!floodplains_glo(:, :) = 0
3891	DO ig = 1,nbp_glo
3892	DO ib = 1,nbasmax
3893	floodplains_glo(ig, ib) = routing_area_glo(ig,ib) * floodplains_glo(ig, ib)
3894	END DO
3895	END DO
3896	END IF
3897	!
3898	! Verifications of the routing graph.
3899	!
3900	nbhtumon = 0
3901	DO ig = 1,nbp_glo
3902	! Noramlize the areas so that differences in precision of area compution by RoutingPP do not affect the model
3903	!
3904	corr = contfrac_g(ig)*area_g(ig)/SUM(routing_area_glo(ig,:))
3905	IF (ABS(1 - corr) > 0.0002 ) THEN
3906	WRITE(,) "Correcting the HTU area to take into account contfrac", corr
3907	IF ( ABS(1 - corr) > 0.1) THEN
3908	WRITE(,) "Coordinates : ", lalo_g(ig,1), lalo_g(ig,2)
3909	WRITE(,) "Contfrac and area in model : ", contfrac_g(ig), area_g(ig)
3910	WRITE(,) "Total grid area in graph file : ", SUM(routing_area_glo(ig,:))
3911	WRITE(,) "The new areas are : ", SUM(routing_area_glo(ig,:)), contfrac_g(ig)*area_g(ig)
3912	WRITE(,) "Correction factor : ", corr
3913	CALL ipslerr(3,'routing_hr_basins_p',&
3914	'There is a mismatch in the area of the grid', &
3915	'Either there are issues with the projection of the grid ',' or contfrac mismatches.')
3916	ELSE
3917	CALL ipslerr(2,'routing_hr_basins_p',&
3918	'The area of the grid had to be adjusted by less than 10% :', &
3919	' ',' ')
3920	ENDIF
3921	ENDIF
3922	DO ib = 1,nbasmax
3923	routing_area_glo(ig,ib) = corr*routing_area_glo(ig,ib)
3924	ENDDO
3925	!
3926	!
3927	DO ib = 1,nbasmax
3928	!
3929	IF (topo_resid_glo(ig,ib) <= zero .AND. route_tobasin_glo(ig, ib) .LE. nbasmax+3) THEN
3930	! If the basin has no surface we change silently as it does not matter.
3931	IF ( routing_area_glo(ig,ib) > zero ) THEN
3932	CALL ipslerr(2,'routing_hr_basins_p',&
3933	'Some zero topo_resid (topoindex) values were encoutered and replaced here :', &
3934	' ',' ')
3935	WRITE(,) "routing_hr_basins_p : topo_resid_glo : ", topo_resid_glo(ig,ib), routing_area_glo(ig,ib)
3936	WRITE(,) "routing_hr_basins_p : Coordinates : ", lalo_g(ig,1), lalo_g(ig,2)
3937	topo_resid_glo(ig,ib) = 10
3938	stream_resid_glo(ig,ib) = 10
3939	WRITE(,) "routing_hr_basins_p : New topo_resid_glo : ", topo_resid_glo(ig,ib)
3940	ELSE
3941	topo_resid_glo(ig,ib) = 10
3942	stream_resid_glo(ig,ib) = 10
3943	ENDIF
3944	ENDIF
3945	!
3946	!
3947	IF ( route_togrid_glo(ig, ib) > nbp_glo ) THEN
3948	IF ( route_tobasin_glo(ig,ib) <= nbasmax+3 ) THEN
3949	WRITE(,) "Issues with the global grid : ", ig, ib, route_togrid_glo(ig, ib), route_tobasin_glo(ig,ib)
3950	CALL ipslerr(3,'routing_hr_basins_p','route_togrid is not compatible with the model configuration', &
3951	' ',' ')
3952	ENDIF
3953	ELSE
3954	ic = 0
3955	ign = ig
3956	ibn = ib
3957	! Locate outflow point
3958	DO WHILE (ibn .GT. 0 .AND. ibn .LE. nbasmax .AND. ic .LT. nbasmax*nbp_glo)
3959	ic = ic + 1
3960	og = ign
3961	ob = ibn
3962	ign = route_togrid_glo(og, ob)
3963	ibn = route_tobasin_glo(og, ob)
3964	!
3965	IF (ibn .GT. nbasmax+3 .OR. ign .GT. nbp_glo) THEN
3966	WRITE(,) "Reached point ", ign, ibn, " on condition ", nbasmax+3, nbp_glo
3967	WRITE(,) "Why do we flow into basin :", route_tobasin_glo(og, ob), " at ", og,ob
3968	WRITE(,) "Coordinates : ", lalo_g(ob,1), lalo_g(ob,2)
3969	WRITE(,) "neighbours_g : ", MINVAL(neighbours_g(ob,:))
3970	CALL ipslerr(3,'routing_hr_basins_p','The river flows into a place outside of the grid.', &
3971	' ',' ')
3972	ENDIF
3973	ENDDO
3974	IF ( ic .GE. nbasmax*nbp_glo) THEN
3975	WRITE(,) "Some river did not converge on point ", ig, ib, ic
3976	WRITE(,) "The start point in the graph was : ", lalo_g(ig,2), lalo_g(ig,1), ib
3977	WRITE(,) "The last point we passed through was : ", lalo_g(og,2), lalo_g(og,1), ob
3978	WRITE(,) "The next one would be : ", lalo_g(ign,2), lalo_g(ign,1), ibn
3979	og = route_togrid_glo(ign, ibn)
3980	ob = route_tobasin_glo(ign, ibn)
3981	WRITE(,) "The after next HTU would be : ", lalo_g(og,2), lalo_g(og,1), ob
3982	WRITE(,) "Last information : ", ign, ibn
3983	CALL ipslerr(3,'routing_hr_basins_p','The river never flows into an outflow point.', &
3984	' ',' ')
3985	ENDIF
3986	ENDIF
3987	ENDDO
3988	!
3989	! Count stations to be monitored
3990	!
3991	DO im=1,nbasmon
3992	IF ( HTUdiag_glo(ig,im) > 0 ) THEN
3993	nbhtumon = nbhtumon + 1
3994	ENDIF
3995	ENDDO
3996	ENDDO
3997	WRITE(numout,*) "Found a total of ", nbhtumon, " HTUs to be monitored and written into HTUhgmon"
3998	!
3999	! Compute num_largest
4000	!
4001	num_largest = COUNT(route_tobasin_glo .EQ. nbasmax+3)
4002	WRITE(numout,*) "After _basins_p : Number of largest rivers : ", COUNT(route_tobasin_glo .EQ. nbasmax+3)
4003	ENDIF
4004	!
4005	CALL bcast(num_largest)
4006	CALL bcast(nbasmax)
4007	CALL bcast(nbasmon)
4008	CALL bcast(inflows)
4009	!
4010	CALL scatter(routing_area_glo,routing_area_loc)
4011	IF ( do_floodplains ) THEN
4012	CALL scatter(floodplains_glo,floodplains_loc)
4013	CALL scatter(floodcri_glo, floodcri_loc)
4014	CALL scatter(fp_beta_glo, fp_beta_loc)
4015	END IF
4016	CALL scatter(global_basinid_glo, global_basinid_loc)
4017	CALL scatter(topo_resid_glo, topo_resid_loc)
4018	CALL scatter(stream_resid_glo, stream_resid_loc)
4019	CALL scatter(route_togrid_glo, route_togrid_loc)
4020	CALL scatter(route_tobasin_glo, route_tobasin_loc)
4021	CALL scatter(route_nbintobas_glo, route_nbintobas_loc)
4022	CALL scatter(hydrodiag_glo, hydrodiag_loc)
4023	CALL scatter(HTUdiag_glo, HTUdiag_loc)
4024	IF ( do_floodplains .AND. dofloodoverflow ) THEN
4025	CALL scatter(orog_min_glo, orog_min_loc)
4026	END IF
4027	!
4028	CALL bcast(stream_tcst)
4029	CALL bcast(fast_tcst)
4030	CALL bcast(slow_tcst)
4031	CALL bcast(flood_tcst)
4032	CALL bcast(swamp_cst)
4033	CALL bcast(lim_floodcri)
4034	CALL bcast(stream_maxresid)
4035	!
4036	END SUBROUTINE routing_hr_basins_p
4037	!
4038	!! ================================================================================================================================
4039	!! SUBROUTINE : routing_hr_graphinfo
4040	!!
4041	!>\BRIEF Extract some basic information from the routing graph file which cannot be obtained through IOIPSL.
4042	!!
4043	!! ================================================================================================================================
4044	SUBROUTINE routing_hr_graphinfo(filename, basmax, infmax, basmon, undef, tstream, tfast, tslow, tflood, cswamp, lfpcri)
4045	!
4046	USE netcdf
4047	!
4048	IMPLICIT NONE
4049	!
4050	!! 0. Variables and parameter declaration
4051	!! 0.1 Input variables
4052	CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open
4053	INTEGER(i_std), INTENT(inout) :: basmax !! maximum number of HTUs
4054	INTEGER(i_std), INTENT(inout) :: basmon !! Number of HTUs to be monitored by grid box.
4055	INTEGER(i_std), INTENT(inout) :: infmax !! Maximum number of inflows.
4056	REAL(r_std), INTENT(out) :: undef
4057	REAL(r_std), INTENT(out) :: tstream, tfast, tslow, tflood, cswamp !! Time constants to be extracted
4058	REAL(r_std), INTENT(out) :: lfpcri !! Constant lim_floodcri to be taken from graph file.
4059	!
4060	INTEGER(i_std) :: rcode, nid, dimid, ndims, nvars
4061	INTEGER(i_std), DIMENSION(4) :: dimids
4062	INTEGER(i_std) :: iv, ndimsvar
4063	CHARACTER(LEN=20) :: dname, varname
4064	!
4065	!
4066	IF (is_root_prc) THEN
4067	!
4068	rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
4069	IF (rcode == NF90_NOERR) THEN
4070	!
4071	! Get graph file version
4072	!
4073	rcode = nf90_get_att(nid, NF90_GLOBAL, "RoutingPPVersion", graphfile_version)
4074	IF (rcode /= NF90_NOERR) THEN
4075	graphfile_version = 0.0
4076	ENDIF
4077	!
4078	! Assumes that the number of HTUs is in the dimension z
4079	!
4080	rcode = nf90_inq_dimid(nid, "z", dimid)
4081	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inq_dimid for z', &
4082	TRIM(nf90_strerror(rcode)),'')
4083	rcode = nf90_inquire_dimension(nid, dimid, dname, basmax)
4084	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire_dimension basmax', &
4085	TRIM(nf90_strerror(rcode)),'')
4086	!
4087	rcode = nf90_inq_dimid(nid, "inflow", dimid)
4088	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inq_dimid for inflow', &
4089	TRIM(nf90_strerror(rcode)),'')
4090	rcode = nf90_inquire_dimension(nid, dimid, dname, infmax)
4091	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire_dimension inflows', &
4092	TRIM(nf90_strerror(rcode)),'')
4093	!
4094	rcode = nf90_inq_dimid(nid, "htumon", dimid)
4095	IF (rcode /= NF90_NOERR) THEN
4096	MonitoringinGraph = .FALSE.
4097	basmon = 1
4098	ELSE
4099	rcode = nf90_inquire_dimension(nid, dimid, dname, basmon)
4100	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire_dimension for basmon', &
4101	TRIM(nf90_strerror(rcode)),'')
4102	MonitoringinGraph = .TRUE.
4103	ENDIF
4104	!
4105	!
4106	rcode = NF90_INQUIRE (nid, nDimensions=ndims, nVariables=nvars)
4107	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire', &
4108	TRIM(nf90_strerror(rcode)),'')
4109	!
4110	DO iv=1,nvars
4111	!
4112	rcode = NF90_INQUIRE_VARIABLE(nid, iv, name=varname, ndims=ndimsvar, dimids=dimids)
4113	!
4114	SELECT CASE (varname)
4115	CASE ("basin_area")
4116	rcode = NF90_GET_ATT(nid, iv, "missing_value", undef)
4117	IF (rcode /= NF90_NOERR) THEN
4118	IF ( rcode == NF90_ENOTATT ) THEN
4119	rcode = NF90_GET_ATT(nid, iv, "_FillValue", undef)
4120	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Did not get FillValue with nf90_get_att', &
4121	TRIM(nf90_strerror(rcode)),'')
4122	ELSE
4123	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_att', &
4124	TRIM(nf90_strerror(rcode)),'')
4125	ENDIF
4126	ENDIF
4127	CASE("StreamTimeCst")
4128	rcode = NF90_GET_VAR(nid,iv,tstream)
4129	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable StreamTimeCst', '','')
4130	! If in an old version convert 10^3d/km into s/km
4131	IF (graphfile_version < 1.0) THEN
4132	tstream = tstream/1000*one_day
4133	ENDIF
4134	CASE("FastTimeCst")
4135	rcode = NF90_GET_VAR(nid,iv,tfast)
4136	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable FastTimeCst', '','')
4137	! If in an old version convert 10^3d/km into s/km
4138	IF (graphfile_version < 1.0) THEN
4139	tfast = tfast/1000*one_day
4140	ENDIF
4141	CASE("SlowTimeCst")
4142	rcode = NF90_GET_VAR(nid,iv,tslow)
4143	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable SlowTimeCst', '','')
4144	! If in an old version convert 10^3d/km into s/km
4145	IF (graphfile_version < 1.0) THEN
4146	tslow = tslow/1000*one_day
4147	ENDIF
4148	CASE("FloodTimeCst")
4149	rcode = NF90_GET_VAR(nid,iv,tflood)
4150	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable FloodTimeCst', '','')
4151	! If in an old version convert 10^3d/km into s/km
4152	IF (graphfile_version < 1.0) THEN
4153	tflood = tflood/1000*one_day
4154	ENDIF
4155	CASE("SwampCst")
4156	rcode = NF90_GET_VAR (nid,iv,cswamp)
4157	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable SwampCst', '','')
4158	CASE("MaxTimeStep")
4159	rcode = NF90_GET_VAR (nid,iv,maxtimestep)
4160	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable MaxTimeStep', '','')
4161	CASE("LimFloodcri")
4162	rcode = NF90_GET_VAR(nid,iv,lfpcri)
4163	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable LimFloodcri', '','')
4164	END SELECT
4165	ENDDO
4166	rcode = NF90_CLOSE(nid)
4167	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_close', &
4168	TRIM(nf90_strerror(rcode)),'')
4169	!
4170	! Before RoutingGraph version 2.5 the lfpcri parameter was hardcoded in routing.f90 and set to 2m.
4171	! This information is preserved here.
4172	IF ( graphfile_version < 2.5 ) THEN
4173	lfpcri = 2.0
4174	ENDIF
4175	!
4176	ELSE
4177	! Case without Graphfile. So we consider that the information will be found in the restart.
4178	CALL ipslerr_p(2, 'routing_hr_graphinfo', 'Could not open the rotung_graph.nc file', &
4179	"Expect to find all the information needed in the restart file.",'')
4180	!
4181	MonitoringinGraph = .FALSE.
4182	basmax = -1
4183	basmon = -1
4184	infmax = -1
4185	undef = undef_sechiba
4186	tstream = -1
4187	tfast = -1
4188	tslow = -1
4189	tflood = -1
4190	cswamp = -1
4191	lfpcri = -1
4192	ENDIF
4193	ENDIF
4194	!!
4195	CALL bcast(MonitoringinGraph)
4196	CALL bcast(basmax)
4197	CALL bcast(basmon)
4198	CALL bcast(infmax)
4199	CALL bcast(undef)
4200	CALL bcast(tstream)
4201	CALL bcast(tfast)
4202	CALL bcast(tslow)
4203	CALL bcast(tflood)
4204	CALL bcast(cswamp)
4205	CALL bcast(lfpcri)
4206	!!
4207	!!
4208	END SUBROUTINE routing_hr_graphinfo
4209	!
4210	!! ================================================================================================================================
4211	!! SUBROUTINE : routing_hr_indexfilegrid
4212	!!
4213	!>\BRIEF Locates all the points of the routing graph file on the model grid. This ensure that no assumption is made on the
4214	!! orientation of the grid in the file.
4215	!!
4216	!! ================================================================================================================================
4217	SUBROUTINE routing_hr_indexfilegrid(im, jm, nbl, lon, lat, landindex, indextab, il2il)
4218	INTEGER(i_std), INTENT(IN) :: im, jm, nbl
4219	REAL(r_std), INTENT(IN) :: lon(im,jm), lat(im,jm)
4220	REAL(r_std), INTENT(IN) :: landindex(im,jm)
4221	INTEGER(i_std), INTENT(INOUT) :: indextab(im,jm)
4222	INTEGER(i_std), INTENT(OUT) :: il2il(nbl)
4223	!
4224	INTEGER(i_std) :: il,i,j
4225	INTEGER(i_std) :: f(2)
4226	INTEGER(i_std) :: ih, jh, ir, jr
4227	REAL(r_std) :: nd
4228	REAL(r_std), DIMENSION(im,jm) :: dist
4229	REAL(r_std) :: mindist = 1000. !! Minimum distance in m between two points to be matched.
4230	!
4231	indextab(:,:) = -1
4232	ih = INT(im/2.)
4233	ir = NINT(im/2.)+1
4234	jh = INT(jm/2.)
4235	jr = NINT(jm/2.)+1
4236	!
4237	DO il=1,nbl
4238	dist(:,:) = undef_sechiba
4239	DO i=MAX(1,ih-ir),MIN(im,ih+ir)
4240	DO j=MAX(1,jh-jr),MIN(jm,jh+jr)
4241	dist(i,j) = haversine_distance(lon(i,j), lat(i,j), lalo_g(il,2), lalo_g(il,1))
4242	ENDDO
4243	ENDDO
4244	f=MINLOC(dist)
4245	IF ( dist(f(1),f(2)) < mindist ) THEN
4246	indextab(f(1),f(2)) = il
4247	il2il(NINT(landindex(f(1),f(2)))) = il
4248	dist(f(1),f(2)) = undef_sechiba
4249	ELSE
4250	CALL ipslerr(3,'routing_hr_indexfilegrid',&
4251	'Distance of the closest point in the two grids is too large. ', &
4252	'Are you sure the routing graph file is on the correct grid ?',&
4253	' ')
4254	ENDIF
4255	!
4256	! See if the next point is close by
4257	!
4258	nd = haversine_distance(lalo_g(il,2), lalo_g(il,1), lalo_g(MIN(il+1,nbl),2), lalo_g(MIN(il+1,nbl),1))
4259	IF ( nd < MINVAL(dist)*3 ) THEN
4260	! The next point is close so zoom in
4261	ih = f(1)
4262	ir = 4
4263	jh = f(2)
4264	jr = 4
4265	ELSE
4266	! Back to starting conditions as the next point is far away
4267	ih = INT(im/2.)
4268	ir = NINT(im/2.)+1
4269	jh = INT(jm/2.)
4270	jr = NINT(jm/2.)+1
4271	ENDIF
4272	ENDDO
4273	END SUBROUTINE routing_hr_indexfilegrid
4274	!
4275	!! ================================================================================================================================
4276	!! SUBROUTINE : routing_hr_convertlandpts
4277	!!
4278	!>\BRIEF In case the order of land points was different in RoutingPreProc and the model. The route_togrid is corrected.
4279	!!
4280	!! ================================================================================================================================
4281	!
4282	SUBROUTINE routing_hr_convertlandpts(nbl, nbas, land2land, route_togrid)
4283	INTEGER(i_std), INTENT(IN) :: nbl, nbas
4284	INTEGER(i_std), INTENT(IN), DIMENSION(nbl) :: land2land
4285	INTEGER(i_std), INTENT(INOUT), DIMENSION(nbl,nbas) :: route_togrid
4286	!
4287	INTEGER(i_std) :: ip, ib
4288	!
4289	DO ip=1,nbl
4290	DO ib=1,nbas
4291	IF ( route_togrid(ip,ib) < undef_int .AND. route_togrid(ip,ib) > 0 ) THEN
4292	route_togrid(ip,ib) = land2land(route_togrid(ip,ib))
4293	ELSE
4294	route_togrid(ip,ib) = ip
4295	ENDIF
4296	ENDDO
4297	ENDDO
4298	!
4299	END SUBROUTINE routing_hr_convertlandpts
4300	!
4301	!! ================================================================================================================================
4302	!! SUBROUTINE : routing_hr_inflows
4303	!!
4304	!>\BRIEF Calculate the inflows from the outflows information.
4305	!!
4306	!! ================================================================================================================================
4307	!
4308	SUBROUTINE routing_hr_inflows(nbl, nbas, inf, floodplains_glo,route_innum_glo,route_ingrid_glo,route_inbasin_glo)
4309
4310	IMPLICIT None
4311
4312	INTEGER(i_std), INTENT(IN) :: nbl, nbas, inf
4313	REAL(r_std), INTENT(IN), DIMENSION(nbl,nbas) :: floodplains_glo
4314	INTEGER(i_std), INTENT(INOUT), DIMENSION(nbl,nbas) :: route_innum_glo
4315	INTEGER(i_std), INTENT(INOUT), DIMENSION(nbl,nbas, inf) :: route_ingrid_glo, route_inbasin_glo
4316	!
4317	INTEGER(i_std) :: ig, ib, og, ob
4318	!
4319	route_innum_glo(:,:) = 0
4320	route_ingrid_glo(:,:,:) = 0
4321	route_inbasin_glo(:,:,:) = 0
4322	DO ig=1,nbl
4323	DO ib=1,nbas
4324	IF (floodplains_glo(ig,ib) .GT. 0) THEN
4325	og = route_togrid_glo(ig,ib)
4326	ob = route_tobasin_glo(ig,ib)
4327	IF (ob .LE. nbasmax) THEN
4328	IF (floodplains_glo(og,ob) .GT. 0) THEN
4329	route_innum_glo(og, ob) = route_innum_glo(og, ob) + 1
4330	route_ingrid_glo(og,ob,route_innum_glo(og, ob)) = ig
4331	route_inbasin_glo(og,ob,route_innum_glo(og, ob)) = ib
4332	END IF
4333	END IF
4334	END IF
4335	ENDDO
4336	ENDDO
4337	END SUBROUTINE routing_hr_inflows
4338	!
4339	!!
4340	!! ================================================================================================================================
4341	!! SUBROUTINE : routing_hr_landgather
4342	!!
4343	!>\BRIEF Gathers the routing information onto landpoints, i.e. goes from an X/Y grid to a list of land points.
4344	!!
4345	!! ================================================================================================================================
4346	!
4347	SUBROUTINE routing_hr_landgather_r(im,jm,nbas,nbl,indextab,ijfield,landfield,def)
4348	!
4349	INTEGER(i_std), INTENT(IN) :: im,jm,nbas,nbl
4350	INTEGER(i_std), INTENT(IN) :: indextab(im,jm)
4351	REAL(r_std), INTENT(IN), DIMENSION(im,jm,nbas) :: ijfield
4352	REAL(r_std), INTENT(OUT), DIMENSION(nbl,nbas) :: landfield
4353	REAL(r_std), INTENT(IN) :: def
4354	!
4355	INTEGER(i_std) :: i,j,k
4356	!
4357	DO i=1,im
4358	DO j=1,jm
4359	IF ( indextab(i,j) > 0 ) THEN
4360	DO k=1,nbas
4361	! Catch undef or NaN values
4362	IF (ijfield(i,j,k) >= undef_graphfile .OR. ijfield(i,j,k) /= ijfield(i,j,k)) THEN
4363	landfield(indextab(i,j),k) = def
4364	ELSE
4365	landfield(indextab(i,j),k) = ijfield(i,j,k)
4366	ENDIF
4367	ENDDO
4368	ENDIF
4369	ENDDO
4370	ENDDO
4371	END SUBROUTINE routing_hr_landgather_r
4372	SUBROUTINE routing_hr_landgather_i2(im,jm,nbas,nbl,indextab,ijfield,landfield,def)
4373	!
4374	INTEGER(i_std), INTENT(IN) :: im,jm,nbas,nbl
4375	INTEGER(i_std), INTENT(IN) :: indextab(im,jm)
4376	REAL(r_std), INTENT(IN), DIMENSION(im,jm,nbas) :: ijfield
4377	INTEGER(i_std), INTENT(OUT), DIMENSION(nbl,nbas) :: landfield
4378	INTEGER(i_std), INTENT(IN) :: def
4379	!
4380	INTEGER(i_std) :: i,j,in
4381	!
4382	DO i=1,im
4383	DO j=1,jm
4384	IF ( indextab(i,j) > 0 ) THEN
4385	DO in=1,nbas
4386	IF (ijfield(i,j,in) .GE. undef_int ) THEN
4387	landfield(indextab(i,j),in) = def
4388	ELSE
4389	landfield(indextab(i,j),in) = ijfield(i,j,in)
4390	ENDIF
4391	ENDDO
4392	ENDIF
4393	ENDDO
4394	ENDDO
4395	END SUBROUTINE routing_hr_landgather_i2
4396	SUBROUTINE routing_hr_landgather_i1(im,jm,nbl,indextab,ijfield,landfield,def)
4397	!
4398	INTEGER(i_std), INTENT(IN) :: im,jm,nbl
4399	INTEGER(i_std), INTENT(IN) :: indextab(im,jm)
4400	REAL(r_std), INTENT(IN), DIMENSION(im,jm) :: ijfield
4401	INTEGER(i_std), INTENT(OUT), DIMENSION(nbl) :: landfield
4402	INTEGER(i_std), INTENT(IN) :: def
4403	!
4404	INTEGER(i_std) :: i,j
4405	!
4406	DO i=1,im
4407	DO j=1,jm
4408	IF ( indextab(i,j) > 0 ) THEN
4409	IF (ijfield(i,j) .GE. undef_int ) THEN
4410	landfield(indextab(i,j)) = def
4411	ELSE
4412	landfield(indextab(i,j)) = ijfield(i,j)
4413	ENDIF
4414	ENDIF
4415	ENDDO
4416	ENDDO
4417	END SUBROUTINE routing_hr_landgather_i1
4418	!
4419	!
4420	!! ================================================================================================================================
4421	!! SUBROUTINE : routing_hr_irrigmap
4422	!!
4423	!>\BRIEF This subroutine interpolates the 0.5x0.5 degree based map of irrigated areas to the resolution of the model.
4424	!!
4425	!! DESCRIPTION (definitions, functional, design, flags) : None
4426	!!
4427	!! RECENT CHANGE(S): None
4428	!!
4429	!! MAIN OUTPUT VARIABLE(S):
4430	!!
4431	!! REFERENCES : None
4432	!!
4433	!! FLOWCHART : None
4434	!! \n
4435	!_ ================================================================================================================================
4436
4437	SUBROUTINE routing_hr_irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, &
4438	& init_irrig, irrigated, init_flood, init_swamp, swamp, hist_id, hist2_id)
4439	!
4440	IMPLICIT NONE
4441	!
4442	!! PARAMETERS
4443	INTEGER(i_std), PARAMETER :: ilake = 1 !! Number of type of lakes area (unitless)
4444	INTEGER(i_std), PARAMETER :: idam = 2 !! Number of type of dams area (unitless)
4445	INTEGER(i_std), PARAMETER :: iflood = 3 !! Number of type of floodplains area (unitless)
4446	INTEGER(i_std), PARAMETER :: iswamp = 4 !! Number of type of swamps area (unitless)
4447	INTEGER(i_std), PARAMETER :: isal = 5 !! Number of type of salines area (unitless)
4448	INTEGER(i_std), PARAMETER :: ipond = 6 !! Number of type of ponds area (unitless)
4449	INTEGER(i_std), PARAMETER :: ntype = 6 !! Number of types of flooded surfaces (unitless)
4450
4451	!! INPUT VARIABLES
4452	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
4453	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Index on the global map.
4454	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
4455	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
4456	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m)
4457	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
4458	INTEGER(i_std), INTENT(in) :: hist_id !! Access to history file (unitless)
4459	INTEGER(i_std), INTENT(in) :: hist2_id !! Access to history file 2 (unitless)
4460	LOGICAL, INTENT(in) :: init_irrig !! Logical to initialize the irrigation (true/false)
4461	LOGICAL, INTENT(in) :: init_flood !! Logical to initialize the floodplains (true/false)
4462	LOGICAL, INTENT(in) :: init_swamp !! Logical to initialize the swamps (true/false)
4463	!
4464	!! OUTPUT VARIABLES
4465	REAL(r_std), INTENT(out) :: irrigated(:) !! Irrigated surface in each grid box (m^2)
4466	!! REAL(r_std), INTENT(out) :: floodplains(:) !! Surface which can be inundated in each grid box (m^2)
4467	REAL(r_std), INTENT(out) :: swamp(:) !! Surface which can be swamp in each grid box (m^2)
4468	!
4469	!! LOCAL VARIABLES
4470	! Interpolation variables
4471	!
4472	INTEGER(i_std) :: nbpmax, nix, njx, fopt !!
4473	CHARACTER(LEN=30) :: callsign !!
4474	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_lu !! Resolution read on the map
4475	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask !! Mask to exclude some points (unitless)
4476	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrsub_area !! Area on the fine grid (m^2)
4477	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: irrsub_index !! Indices of the points we need on the fine grid (unitless)
4478	INTEGER :: ALLOC_ERR !!
4479	LOGICAL :: ok_interpol = .FALSE. !! Flag for interpolation (true/false)
4480	!
4481	CHARACTER(LEN=80) :: filename !! Name of the netcdf file (unitless)
4482	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, itype !! Indices (unitless)
4483	REAL(r_std) :: lev(1), date, dt, coslat !!
4484	INTEGER(i_std) :: itau(1) !!
4485	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: latrel !! Latitude
4486	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lonrel !! Longitude
4487	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrigated_frac !! Irrigated fraction of the grid box (unitless;0-1)
4488	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: flood_fracmax !! Maximal flooded fraction of the grid box (unitless;0-1)
4489	REAL(r_std) :: area_irrig !! Irrigated surface in the grid box (m^2)
4490	REAL(r_std) :: area_flood(ntype) !! Flooded surface in the grid box (m^2)
4491	!!$ REAL(r_std) :: irrigmap(nbpt)
4492	!!$ REAL(r_std) :: swampmap(nbpt)
4493
4494	!_ ================================================================================================================================
4495
4496	!
4497	!Config Key = IRRIGATION_FILE
4498	!Config Desc = Name of file which contains the map of irrigated areas
4499	!Config Def = floodplains.nc
4500	!Config If = DO_IRRIGATION OR DO_FLOODPLAINS
4501	!Config Help = The name of the file to be opened to read the field
4502	!Config with the area in m^2 of the area irrigated within each
4503	!Config 0.5 0.5 deg grid box. The map currently used is the one
4504	!Config developed by the Center for Environmental Systems Research
4505	!Config in Kassel (1995).
4506	!Config Units = [FILE]
4507	!
4508	filename = 'floodplains.nc'
4509	CALL getin_p('IRRIGATION_FILE',filename)
4510	!
4511	IF (is_root_prc) THEN
4512	CALL flininfo(filename,iml, jml, lml, tml, fid)
4513	CALL flinclo(fid)
4514	ELSE
4515	iml = 0
4516	jml = 0
4517	lml = 0
4518	tml = 0
4519	ENDIF
4520	!
4521	CALL bcast(iml)
4522	CALL bcast(jml)
4523	CALL bcast(lml)
4524	CALL bcast(tml)
4525	!
4526	!
4527	!
4528	ALLOCATE (latrel(iml,jml), STAT=ALLOC_ERR)
4529	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for latrel','','')
4530
4531	ALLOCATE (lonrel(iml,jml), STAT=ALLOC_ERR)
4532	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for lonrel','','')
4533
4534	ALLOCATE (irrigated_frac(iml,jml), STAT=ALLOC_ERR)
4535	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for irrigated_frac','','')
4536
4537	ALLOCATE (flood_fracmax(iml,jml,ntype), STAT=ALLOC_ERR)
4538	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for flood_fracmax','','')
4539
4540	IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lonrel, latrel, lev, tml, itau, date, dt, fid)
4541
4542	CALL bcast(lonrel)
4543	CALL bcast(latrel)
4544	!
4545	IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 0, 0, irrigated_frac)
4546	CALL bcast(irrigated_frac)
4547	IF (is_root_prc) CALL flinget(fid, 'lake', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,ilake))
4548	IF (is_root_prc) CALL flinget(fid, 'dam', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,idam))
4549	IF (is_root_prc) CALL flinget(fid, 'flood', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,iflood))
4550	IF (is_root_prc) CALL flinget(fid, 'swamp', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,iswamp))
4551	IF (is_root_prc) CALL flinget(fid, 'saline', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,isal))
4552	IF (is_root_prc) CALL flinget(fid, 'pond', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,ipond))
4553	CALL bcast(flood_fracmax)
4554	!
4555	IF (is_root_prc) CALL flinclo(fid)
4556	!
4557	! Set to zero all fraction which are less than 0.5%
4558	!
4559	DO ip=1,iml
4560	DO jp=1,jml
4561	!
4562	IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN
4563	irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100.
4564	IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero
4565	ENDIF
4566	!
4567	DO itype=1,ntype
4568	IF ( flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN
4569	flood_fracmax(ip,jp,itype) = flood_fracmax(ip,jp,itype)/100
4570	IF ( flood_fracmax(ip,jp,itype) < 0.005 ) flood_fracmax(ip,jp,itype) = zero
4571	ENDIF
4572	ENDDO
4573	!
4574	ENDDO
4575	ENDDO
4576
4577	IF (printlev>=2) THEN
4578	WRITE(numout,*) 'lonrel : ', MAXVAL(lonrel), MINVAL(lonrel)
4579	WRITE(numout,*) 'latrel : ', MAXVAL(latrel), MINVAL(latrel)
4580	WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), &
4581	MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba)
4582	WRITE(numout,*) 'flood_fracmax : ', MINVAL(flood_fracmax, MASK=flood_fracmax .GT. 0), &
4583	MAXVAL(flood_fracmax, MASK=flood_fracmax .LT. undef_sechiba)
4584	END IF
4585
4586	! Consider all points a priori
4587	!
4588	ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR)
4589	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for resol_lu','','')
4590
4591	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4592	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for mask','','')
4593	mask(:,:) = 0
4594
4595	DO ip=1,iml
4596	DO jp=1,jml
4597	!
4598	! Exclude the points where we are close to the missing value.
4599	!
4600	!MG This condition cannot be applied in floodplains/swamps configuration because
4601	! the same mask would be used for the interpolation of irrigation, floodplains and swamps maps.
4602	! IF ( irrigated_frac(ip,jp) < undef_sechiba ) THEN
4603	mask(ip,jp) = 1
4604	! ENDIF
4605	!
4606	! Resolution in longitude
4607	!
4608	coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos )
4609	IF ( ip .EQ. 1 ) THEN
4610	resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip,jp) ) * pi/180. * R_Earth * coslat
4611	ELSEIF ( ip .EQ. iml ) THEN
4612	resol_lu(ip,jp,1) = ABS( lonrel(ip,jp) - lonrel(ip-1,jp) ) * pi/180. * R_Earth * coslat
4613	ELSE
4614	resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip-1,jp) )/2. * pi/180. * R_Earth * coslat
4615	ENDIF
4616	!
4617	! Resolution in latitude
4618	!
4619	IF ( jp .EQ. 1 ) THEN
4620	resol_lu(ip,jp,2) = ABS( latrel(ip,jp) - latrel(ip,jp+1) ) * pi/180. * R_Earth
4621	ELSEIF ( jp .EQ. jml ) THEN
4622	resol_lu(ip,jp,2) = ABS( latrel(ip,jp-1) - latrel(ip,jp) ) * pi/180. * R_Earth
4623	ELSE
4624	resol_lu(ip,jp,2) = ABS( latrel(ip,jp-1) - latrel(ip,jp+1) )/2. * pi/180. * R_Earth
4625	ENDIF
4626	!
4627	ENDDO
4628	ENDDO
4629	!
4630	! The number of maximum vegetation map points in the GCM grid is estimated.
4631	! Some lmargin is taken.
4632	!
4633	callsign = 'Irrigation map'
4634	ok_interpol = .FALSE.
4635	IF (is_root_prc) THEN
4636	nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2
4637	njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2
4638	nbpmax = nixnjx2
4639	IF (printlev>=1) THEN
4640	WRITE(numout,*) "Projection arrays for ",callsign," : "
4641	WRITE(numout,*) "nbpmax = ",nbpmax, nix, njx
4642	END IF
4643	ENDIF
4644	CALL bcast(nbpmax)
4645
4646	ALLOCATE(irrsub_index(nbpt, nbpmax, 2), STAT=ALLOC_ERR)
4647	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for irrsub_index','','')
4648	irrsub_index(:,:,:)=0
4649
4650	ALLOCATE(irrsub_area(nbpt, nbpmax), STAT=ALLOC_ERR)
4651	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for irrsub_area','','')
4652	irrsub_area(:,:)=zero
4653
4654	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
4655	& iml, jml, lonrel, latrel, mask, callsign, &
4656	& nbpmax, irrsub_index, irrsub_area, ok_interpol)
4657	!
4658	!
4659	WHERE (irrsub_area < 0) irrsub_area=zero
4660	!
4661	! Test here if not all sub_area are larger than 0 if so, then we need to increase nbpmax
4662	!
4663	DO ib=1,nbpt
4664	!
4665	area_irrig = 0.0
4666	area_flood = 0.0
4667	!
4668	DO fopt=1,COUNT(irrsub_area(ib,:) > zero)
4669	!
4670	ip = irrsub_index(ib, fopt, 1)
4671	jp = irrsub_index(ib, fopt, 2)
4672	!
4673	IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN
4674	area_irrig = area_irrig + irrsub_area(ib,fopt)*irrigated_frac(ip,jp)
4675	ENDIF
4676	!
4677	DO itype=1,ntype
4678	IF (flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN
4679	area_flood(itype) = area_flood(itype) + irrsub_area(ib,fopt)*flood_fracmax(ip,jp,itype)
4680	ENDIF
4681	ENDDO
4682	ENDDO
4683	!
4684	! Put the total irrigated and flooded areas in the output variables
4685	!
4686	IF ( init_irrig ) THEN
4687	irrigated(ib) = MIN(area_irrig, resolution(ib,1)resolution(ib,2)contfrac(ib))
4688	IF ( irrigated(ib) < 0 ) THEN
4689	WRITE(numout,*) 'We have a problem here : ', irrigated(ib)
4690	WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
4691	WRITE(numout,*) area_irrig
4692	CALL ipslerr_p(3,'routing_hr_irrigmap','Problem with irrigated...','','')
4693	ENDIF
4694	!!$ ! Compute a diagnostic of the map.
4695	!!$ IF(contfrac(ib).GT.zero) THEN
4696	!!$ irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)resolution(ib,2)contfrac(ib) )
4697	!!$ ELSE
4698	!!$ irrigmap (ib) = zero
4699	!!$ ENDIF
4700	!
4701	ENDIF
4702	!
4703	!
4704	!
4705	IF ( init_swamp ) THEN
4706	swamp(ib) = MIN(area_flood(iswamp), resolution(ib,1)resolution(ib,2)contfrac(ib))
4707	IF ( swamp(ib) < 0 ) THEN
4708	WRITE(numout,*) 'We have a problem here : ', swamp(ib)
4709	WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
4710	WRITE(numout,*) area_flood
4711	CALL ipslerr_p(3,'routing_hr_irrigmap','Problem with swamp...','','')
4712	ENDIF
4713	!!$ ! Compute a diagnostic of the map.
4714	!!$ IF(contfrac(ib).GT.zero) THEN
4715	!!$ swampmap(ib) = swamp(ib) / ( resolution(ib,1)resolution(ib,2)contfrac(ib) )
4716	!!$ ELSE
4717	!!$ swampmap(ib) = zero
4718	!!$ ENDIF
4719	ENDIF
4720	!
4721	!
4722	ENDDO
4723	!
4724	!
4725
4726	IF (printlev>=1) THEN
4727	IF ( init_irrig ) WRITE(numout,*) "Diagnostics irrigated :", MINVAL(irrigated), MAXVAL(irrigated)
4728	IF ( init_flood ) WRITE(numout,*) "Diagnostics floodplains :", MINVAL(floodplains), MAXVAL(floodplains)
4729	IF ( init_swamp ) WRITE(numout,*) "Diagnostics swamp :", MINVAL(swamp), MAXVAL(swamp)
4730	END IF
4731
4732	! No compensation is done for overlapping floodplains, swamp and irrig. At least overlapping will not
4733	! happen between floodplains and swamp alone
4734	! IF ( init_irrig .AND. init_flood ) THEN
4735	! DO ib = 1, nbpt
4736	! surp = (floodplains(ib)+swamp(ib)+irrigated(ib)) / (resolution(ib,1)resolution(ib,2)contfrac(ib))
4737	! IF ( surp .GT. un ) THEN
4738	! floodplains(ib) = floodplains(ib) / surp
4739	! swamp(ib) = swamp(ib) / surp
4740	! irrigated(ib) = irrigated(ib) / surp
4741	! ENDIF
4742	! ENDDO
4743	! ENDIF
4744	!
4745	DEALLOCATE (irrsub_area)
4746	DEALLOCATE (irrsub_index)
4747	!
4748	DEALLOCATE (mask)
4749	DEALLOCATE (resol_lu)
4750	!
4751	DEALLOCATE (lonrel)
4752	DEALLOCATE (latrel)
4753	!
4754	END SUBROUTINE routing_hr_irrigmap
4755	!
4756	!! ================================================================================================================================
4757	!! SUBROUTINE : routing_hr_waterbal
4758	!!
4759	!>\BRIEF This subroutine checks the water balance in the routing module.
4760	!!
4761	!! DESCRIPTION (definitions, functional, design, flags) : None
4762	!!
4763	!! RECENT CHANGE(S): None
4764	!!
4765	!! MAIN OUTPUT VARIABLE(S):
4766	!!
4767	!! REFERENCES : None
4768	!!
4769	!! FLOWCHART : None
4770	!! \n
4771	!_ ================================================================================================================================
4772
4773	SUBROUTINE routing_hr_waterbal(nbpt, reinit, floodout, runoff, drainage, returnflow, &
4774	& reinfiltration, irrigation, riverflow, coastalflow)
4775	!
4776	IMPLICIT NONE
4777	!
4778	!! INPUT VARIABLES
4779	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
4780	LOGICAL, INTENT(in) :: reinit !! Controls behaviour (true/false)
4781	REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt)
4782	REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff (kg/m^2/dt)
4783	REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt)
4784	REAL(r_std), INTENT(in) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box.
4785	!! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt)
4786	REAL(r_std), INTENT(in) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
4787	REAL(r_std), INTENT(in) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
4788	REAL(r_std), INTENT(in) :: riverflow(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
4789	REAL(r_std), INTENT(in) :: coastalflow(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
4790	!
4791	! We sum-up all the water we have in the warious reservoirs
4792	!
4793	REAL(r_std), SAVE :: totw_flood !! Sum of all the water amount in the floodplains reservoirs (kg)
4794	!$OMP THREADPRIVATE(totw_flood)
4795	REAL(r_std), SAVE :: totw_stream !! Sum of all the water amount in the stream reservoirs (kg)
4796	!$OMP THREADPRIVATE(totw_stream)
4797	REAL(r_std), SAVE :: totw_fast !! Sum of all the water amount in the fast reservoirs (kg)
4798	!$OMP THREADPRIVATE(totw_fast)
4799	REAL(r_std), SAVE :: totw_slow !! Sum of all the water amount in the slow reservoirs (kg)
4800	!$OMP THREADPRIVATE(totw_slow)
4801	REAL(r_std), SAVE :: totw_lake !! Sum of all the water amount in the lake reservoirs (kg)
4802	!$OMP THREADPRIVATE(totw_lake)
4803	REAL(r_std), SAVE :: totw_pond !! Sum of all the water amount in the pond reservoirs (kg)
4804	!$OMP THREADPRIVATE(totw_pond)
4805	REAL(r_std), SAVE :: totw_in !! Sum of the water flow in to the routing scheme
4806	!$OMP THREADPRIVATE(totw_in)
4807	REAL(r_std), SAVE :: totw_out !! Sum of the water flow out to the routing scheme
4808	!$OMP THREADPRIVATE(totw_out)
4809	REAL(r_std), SAVE :: totw_return !!
4810	!$OMP THREADPRIVATE(totw_return)
4811	REAL(r_std), SAVE :: totw_irrig !!
4812	!$OMP THREADPRIVATE(totw_irrig)
4813	REAL(r_std), SAVE :: totw_river !!
4814	!$OMP THREADPRIVATE(totw_river)
4815	REAL(r_std), SAVE :: totw_coastal !!
4816	!$OMP THREADPRIVATE(totw_coastal)
4817	REAL(r_std) :: totarea !! Total area of basin (m^2)
4818	REAL(r_std) :: area !! Total area of routing (m^2)
4819	INTEGER(i_std) :: ig !!
4820	!
4821	! Just to make sure we do not get too large numbers !
4822	!
4823	!! PARAMETERS
4824	REAL(r_std), PARAMETER :: scaling = 1.0E+6 !!
4825	REAL(r_std), PARAMETER :: allowed_err = 50. !!
4826
4827	!_ ================================================================================================================================
4828	!
4829	IF ( reinit ) THEN
4830	!
4831	totw_flood = zero
4832	totw_stream = zero
4833	totw_fast = zero
4834	totw_slow = zero
4835	totw_lake = zero
4836	totw_pond = zero
4837	totw_in = zero
4838	!
4839	DO ig=1,nbpt
4840	!
4841	totarea = SUM(routing_area(ig,:))
4842	!
4843	totw_flood = totw_flood + SUM(flood_reservoir(ig,:)/scaling)
4844	totw_stream = totw_stream + SUM(stream_reservoir(ig,:)/scaling)
4845	totw_fast = totw_fast + SUM(fast_reservoir(ig,:)/scaling)
4846	totw_slow = totw_slow + SUM(slow_reservoir(ig,:)/scaling)
4847	totw_lake = totw_lake + lake_reservoir(ig)/scaling
4848	totw_pond = totw_pond + pond_reservoir(ig)/scaling
4849	!
4850	totw_in = totw_in + (runoff(ig)totarea + drainage(ig)totarea - floodout(ig)*totarea)/scaling
4851	!
4852	ENDDO
4853	!
4854	ELSE
4855	!
4856	totw_out = zero
4857	totw_return = zero
4858	totw_irrig = zero
4859	totw_river = zero
4860	totw_coastal = zero
4861	area = zero
4862	!
4863	DO ig=1,nbpt
4864	!
4865	totarea = SUM(routing_area(ig,:))
4866	!
4867	totw_flood = totw_flood - SUM(flood_reservoir(ig,:)/scaling)
4868	totw_stream = totw_stream - SUM(stream_reservoir(ig,:)/scaling)
4869	totw_fast = totw_fast - SUM(fast_reservoir(ig,:)/scaling)
4870	totw_slow = totw_slow - SUM(slow_reservoir(ig,:)/scaling)
4871	totw_lake = totw_lake - lake_reservoir(ig)/scaling
4872	totw_pond = totw_pond - pond_reservoir(ig)/scaling
4873	!
4874	totw_return = totw_return + (reinfiltration(ig)+returnflow(ig))*totarea/scaling
4875	totw_irrig = totw_irrig + irrigation(ig)*totarea/scaling
4876	totw_river = totw_river + riverflow(ig)/scaling
4877	totw_coastal = totw_coastal + coastalflow(ig)/scaling
4878	!
4879	area = area + totarea
4880	!
4881	ENDDO
4882	totw_out = totw_return + totw_irrig + totw_river + totw_coastal
4883	!
4884	! Now we have all the information to balance our water
4885	!
4886	IF ( ABS((totw_flood + totw_stream + totw_fast + totw_slow + totw_lake + totw_pond) - &
4887	& (totw_out - totw_in)) > allowed_err ) THEN
4888	WRITE(numout,*) 'WARNING : Water not conserved in routing. Limit at ', allowed_err, ' 10^6 kg'
4889	WRITE(numout,*) '--Water-- change : flood stream fast ', totw_flood, totw_stream, totw_fast
4890	WRITE(numout,*) '--Water-- change : slow, lake ', totw_slow, totw_lake
4891	WRITE(numout,*) '--Water>>> change in the routing res. : ', totw_flood + totw_stream + totw_fast + totw_slow + totw_lake
4892	WRITE(numout,*) '--Water input : ', totw_in
4893	WRITE(numout,*) '--Water output : ', totw_out
4894	WRITE(numout,*) '--Water output : return, irrig ', totw_return, totw_irrig
4895	WRITE(numout,*) '--Water output : river, coastal ',totw_river, totw_coastal
4896	WRITE(numout,*) '--Water>>> change by fluxes : ', totw_out - totw_in, ' Diff [mm/dt]: ', &
4897	& ((totw_flood + totw_stream + totw_fast + totw_slow + totw_lake) - (totw_out - totw_in))/area
4898
4899	! Stop the model
4900	CALL ipslerr_p(3, 'routing_hr_waterbal', 'Water is not conserved in routing.','','')
4901	ENDIF
4902	!
4903	ENDIF
4904	!
4905	END SUBROUTINE routing_hr_waterbal
4906	!
4907	!
4908	END MODULE routing_highres

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

Download in other formats: