Context Navigation

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

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

Integrate routing scheme "highres" developped by Jan Polcher, tested and integrated in ORCHIDEE_2_2 by Lucia Rinchiuso. This corresponds to the revision 7574 of perso/lucia.rinchiuso/myORCHIDEE_2_2_r7481. The developements of the routing scheme from branches/ORCHIDEE-ROUTING at revision 7545 are taken into account.

See also ticket #842

Property svn:executable set to *

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

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

Download in other formats: