Context Navigation

source: branches/ORCHIDEE_Quest/ORCHIDEE/src_sechiba/routing_simple.f90 @ 7406

Last change on this file since 7406 was 6102, checked in by josefine.ghattas, 6 years ago

Added routing_simple module and routing_wrapper. No changes in defalut set up.

See ticket #581

Property svn:keywords set to Date Revision HeadURL

File size: 132.9 KB

Line
1	! =================================================================================================================================
2	! MODULE : routing_simple
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: The subroutines in this subroutine is only called when ROUTING_METHOD=simple is set in run.def.
13	!! The method can be used for regular latitude-longitude grid or for unstructured grid.
14	!!
15	!! RECENT CHANGE(S): None
16	!!
17	!! REFERENCE(S) :
18	!!
19	!! SVN :
20	!! $HeadURL$
21	!! $Date$
22	!! $Revision$
23	!! \n
24	!_ ================================================================================================================================
25
26
27
28	! Histoire Salee
29	!---------------
30	! La douce riviere
31	! Sortant de son lit
32	! S'est jetee ma chere
33	! dans les bras mais oui
34	! du beau fleuve
35	!
36	! L'eau coule sous les ponts
37	! Et puis les flots s'emeuvent
38	! - N'etes vous pas au courant ?
39	! Il parait que la riviere
40	! Va devenir mer
41	! Roland Bacri
42	!
43
44
45	MODULE routing_simple
46
47	USE ioipsl
48	USE xios_orchidee
49	USE ioipsl_para
50	USE constantes
51	USE constantes_soil
52	USE pft_parameters
53	USE sechiba_io_p
54	USE interpol_help
55	USE grid
56	USE mod_orchidee_para
57
58
59	IMPLICIT NONE
60	PRIVATE
61	PUBLIC :: routing_simple_main, routing_simple_xios_initialize
62	PUBLIC :: routing_simple_initialize, routing_simple_finalize, routing_simple_clear
63
64	!! PARAMETERS
65	INTEGER(i_std), PARAMETER :: nbasmax=5 !! The maximum number of basins we wish to have per grid box (truncation of the model) (unitless)
66	INTEGER(i_std), SAVE :: nbvmax !! The maximum number of basins we can handle at any time during the generation of the maps (unitless)
67	!$OMP THREADPRIVATE(nbvmax)
68	REAL(r_std), SAVE :: fast_tcst = 3.0 !! Property of the fast reservoir (day/m)
69	!$OMP THREADPRIVATE(fast_tcst)
70	REAL(r_std), SAVE :: slow_tcst = 25.0 !! Property of the slow reservoir (day/m)
71	!$OMP THREADPRIVATE(slow_tcst)
72	REAL(r_std), SAVE :: stream_tcst = 0.24 !! Property of the stream reservoir (day/m)
73	!$OMP THREADPRIVATE(stream_tcst)
74	REAL(r_std), SAVE :: flood_tcst = 4.0 !! Property of the floodplains reservoir (day/m)
75	!$OMP THREADPRIVATE(flood_tcst)
76	REAL(r_std), SAVE :: swamp_cst = 0.2 !! Fraction of the river transport that flows to the swamps (unitless;0-1)
77	!$OMP THREADPRIVATE(swamp_cst)
78	!
79	! Relation between volume and fraction of floodplains
80	!
81	REAL(r_std), SAVE :: beta = 2.0 !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless)
82	!$OMP THREADPRIVATE(beta)
83	REAL(r_std), SAVE :: betap = 0.5 !! Ratio of the basin surface intercepted by ponds and the maximum surface of ponds (unitless;0-1)
84	!$OMP THREADPRIVATE(betap)
85	REAL(r_std), SAVE :: floodcri = 2000.0 !! Potential height for which all the basin is flooded (mm)
86	!$OMP THREADPRIVATE(floodcri)
87	!
88	! Relation between maximum surface of ponds and basin surface, and drainage (mm/j) to the slow_res
89	!
90	REAL(r_std), PARAMETER :: pond_bas = 50.0 !! [DISPENSABLE] - not used
91	REAL(r_std), SAVE :: pondcri = 2000.0 !! Potential height for which all the basin is a pond (mm)
92	!$OMP THREADPRIVATE(pondcri)
93
94	REAL(r_std), PARAMETER :: maxevap_lake = 7.5/86400. !! Maximum evaporation rate from lakes (kg/m^2/s)
95	REAL(r_std),SAVE :: dt_routing !! Routing time step (s)
96	!$OMP THREADPRIVATE(dt_routing)
97	INTEGER(i_std), SAVE :: diagunit = 87 !! Diagnostic file unit (unitless)
98	!$OMP THREADPRIVATE(diagunit)
99
100	! Logicals to control model configuration
101	!
102	LOGICAL, SAVE :: dofloodinfilt = .FALSE. !! Logical to choose if floodplains infiltration is activated or not (true/false)
103	!$OMP THREADPRIVATE(dofloodinfilt)
104	LOGICAL, SAVE :: doswamps = .FALSE. !! Logical to choose if swamps are activated or not (true/false)
105	!$OMP THREADPRIVATE(doswamps)
106	LOGICAL, SAVE :: doponds = .FALSE. !! Logical to choose if ponds are activated or not (true/false)
107	!$OMP THREADPRIVATE(doponds)
108
109	INTEGER(i_std), SAVE :: num_largest !! Number of largest river basins which should be treated as independently as rivers
110	!! (not flow into ocean as diffusion coastal flow) (unitless)
111	!$OMP THREADPRIVATE(num_largest)
112	REAL(r_std), SAVE :: time_counter !! Time counter (s)
113	!$OMP THREADPRIVATE(time_counter)
114	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_loc !! Surface of basin (m^2)
115	!$OMP THREADPRIVATE(routing_area_loc)
116	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_loc !! Topographic index of the retention time (m)
117	!$OMP THREADPRIVATE(topo_resid_loc)
118	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_loc !! Grid into which the basin flows (unitless)
119	!$OMP THREADPRIVATE(route_togrid_loc)
120	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_loc !! Basin in to which the water goes (unitless)
121	!$OMP THREADPRIVATE(route_tobasin_loc)
122	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_loc !! Number of basin into current one (unitless)
123	!$OMP THREADPRIVATE(route_nbintobas_loc)
124	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_loc !! ID of basin (unitless)
125	!$OMP THREADPRIVATE(global_basinid_loc)
126	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_loc !! Variable to diagnose the hydrographs
127	!$OMP THREADPRIVATE(hydrodiag_loc)
128	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydroupbasin_loc !! The area upstream of the gauging station (m^2)
129	!$OMP THREADPRIVATE(hydroupbasin_loc)
130
131	! parallelism
132	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
133	!$OMP THREADPRIVATE(routing_area_glo)
134	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m)
135	!$OMP THREADPRIVATE(topo_resid_glo)
136	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless)
137	!$OMP THREADPRIVATE(route_togrid_glo)
138	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless)
139	!$OMP THREADPRIVATE(route_tobasin_glo)
140	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_glo !! Number of basin into current one (unitless)
141	!$OMP THREADPRIVATE(route_nbintobas_glo)
142	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless)
143	!$OMP THREADPRIVATE(global_basinid_glo)
144	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_glo !! Variable to diagnose the hydrographs
145	!$OMP THREADPRIVATE(hydrodiag_glo)
146	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydroupbasin_glo !! The area upstream of the gauging station (m^2)
147	!$OMP THREADPRIVATE(hydroupbasin_glo)
148
149	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: routing_area !! Surface of basin (m^2)
150	!$OMP THREADPRIVATE(routing_area)
151	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: topo_resid !! Topographic index of the retention time (m)
152	!$OMP THREADPRIVATE(topo_resid)
153	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_togrid !! Grid into which the basin flows (unitless)
154	!$OMP THREADPRIVATE(route_togrid)
155	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_tobasin !! Basin in to which the water goes (unitless)
156	!$OMP THREADPRIVATE(route_tobasin)
157	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_nbintobas !! Number of basin into current one (unitless)
158	!$OMP THREADPRIVATE(route_nbintobas)
159	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: global_basinid !! ID of basin (unitless)
160	!$OMP THREADPRIVATE(global_basinid)
161	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: hydrodiag !! Variable to diagnose the hydrographs
162	!$OMP THREADPRIVATE(hydrodiag)
163	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slowflow_diag !! Diagnostic slow flow hydrographs (kg/dt)
164	!$OMP THREADPRIVATE(slowflow_diag)
165	REAL(r_std), SAVE, POINTER, DIMENSION(:) :: hydroupbasin !! The area upstream of the gauging station (m^2)
166	!$OMP THREADPRIVATE(hydroupbasin)
167
168	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigated !! Area equipped for irrigation in each grid box (m^2)
169	!$OMP THREADPRIVATE(irrigated)
170	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodplains !! Maximal surface which can be inundated in each grid box (m^2)
171	!$OMP THREADPRIVATE(floodplains)
172	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: swamp !! Maximal surface of swamps in each grid box (m^2)
173	!$OMP THREADPRIVATE(swamp)
174	!
175	! The reservoirs, variables kept in the restart file.
176	!
177	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_reservoir !! Water amount in the fast reservoir (kg)
178	!$OMP THREADPRIVATE(fast_reservoir)
179	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: slow_reservoir !! Water amount in the slow reservoir (kg)
180	!$OMP THREADPRIVATE(slow_reservoir)
181	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: stream_reservoir !! Water amount in the stream reservoir (kg)
182	!$OMP THREADPRIVATE(stream_reservoir)
183	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir !! Water amount in the floodplains reservoir (kg)
184	!$OMP THREADPRIVATE(flood_reservoir)
185	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_reservoir !! Water amount in the lake reservoir (kg)
186	!$OMP THREADPRIVATE(lake_reservoir)
187	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_reservoir !! Water amount in the pond reservoir (kg)
188	!$OMP THREADPRIVATE(pond_reservoir)
189	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas !! Flooded fraction per basin (unitless;0-1)
190	!$OMP THREADPRIVATE(flood_frac_bas)
191	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_frac !! Pond fraction per grid box (unitless;0-1)
192	!$OMP THREADPRIVATE(pond_frac)
193	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_height !! Floodplain height (mm)
194	!$OMP THREADPRIVATE(flood_height)
195	!
196	! The accumulated fluxes.
197	!
198	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodout_mean !! Accumulated flow out of floodplains (kg/m^2/dt)
199	!$OMP THREADPRIVATE(floodout_mean)
200	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: runoff_mean !! Accumulated runoff (kg/m^2/dt)
201	!$OMP THREADPRIVATE(runoff_mean)
202	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: drainage_mean !! Accumulated drainage (kg/m^2/dt)
203	!$OMP THREADPRIVATE(drainage_mean)
204	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: transpot_mean !! Mean potential transpiration from the plants (kg/m^2/dt)
205	!$OMP THREADPRIVATE(transpot_mean)
206	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: precip_mean !! Accumulated precipitation (kg/m^2/dt)
207	!$OMP THREADPRIVATE(precip_mean)
208	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: humrel_mean !! Mean soil moisture stress, mean root extraction potential (unitless)
209	!$OMP THREADPRIVATE(humrel_mean)
210	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: totnobio_mean !! Mean last total fraction of no bio (unitless;0-1)
211	!$OMP THREADPRIVATE(totnobio_mean)
212	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: vegtot_mean !! Mean potentially vegetated fraction (unitless;0-1)
213	!$OMP THREADPRIVATE(vegtot_mean)
214	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: k_litt_mean !! Mean averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
215	!$OMP THREADPRIVATE(k_litt_mean)
216	!
217	! The averaged outflow fluxes.
218	!
219	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lakeinflow_mean !! Mean lake inflow (kg/m^2/dt)
220	!$OMP THREADPRIVATE(lakeinflow_mean)
221	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: returnflow_mean !! Mean water flow from lakes and swamps which returns to the grid box.
222	!! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
223	!$OMP THREADPRIVATE(returnflow_mean)
224	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: reinfiltration_mean !! Mean water flow which returns to the grid box (kg/m^2/dt)
225	!$OMP THREADPRIVATE(reinfiltration_mean)
226	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigation_mean !! Mean irrigation flux.
227	!! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
228	!$OMP THREADPRIVATE(irrigation_mean)
229	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean !! Mean Outflow of the major rivers.
230	!! The flux will be located on the continental grid but this should be a coastal point (kg/dt)
231	!$OMP THREADPRIVATE(riverflow_mean)
232	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: coastalflow_mean !! Mean outflow on coastal points by small basins.
233	!! This is the water which flows in a disperse way into the ocean (kg/dt)
234	!$OMP THREADPRIVATE(coastalflow_mean)
235	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodtemp !! Temperature to decide if floodplains work (K)
236	!$OMP THREADPRIVATE(floodtemp)
237	INTEGER(i_std), SAVE :: floodtemp_lev !! Temperature level to decide if floodplains work (K)
238	!$OMP THREADPRIVATE(floodtemp_lev)
239	!
240	! Diagnostic variables ... well sort of !
241	!
242	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_netereq !! Irrigation requirement (water requirements by the crop for its optimal growth (kg/m^2/dt)
243	!$OMP THREADPRIVATE(irrig_netereq)
244	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrographs !! Hydrographs at the outflow of the grid box for major basins (kg/dt)
245	!$OMP THREADPRIVATE(hydrographs)
246	!
247	! Diagnostics for the various reservoirs we use (Kg/m^2)
248	!
249	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: fast_diag !! Diagnostic for the fast reservoir (kg/m^2)
250	!$OMP THREADPRIVATE(fast_diag)
251	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slow_diag !! Diagnostic for the slow reservoir (kg/m^2)
252	!$OMP THREADPRIVATE(slow_diag)
253	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: stream_diag !! Diagnostic for the stream reservoir (kg/m^2)
254	!$OMP THREADPRIVATE(stream_diag)
255	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_diag !! Diagnostic for the floodplain reservoir (kg/m^2)
256	!$OMP THREADPRIVATE(flood_diag)
257	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_diag !! Diagnostic for the pond reservoir (kg/m^2)
258	!$OMP THREADPRIVATE(pond_diag)
259	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_diag !! Diagnostic for the lake reservoir (kg/m^2)
260	!$OMP THREADPRIVATE(lake_diag)
261	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: delsurfstor !! Diagnostic of the change in surface water storage (flood, pond and lake reservoirs) (kg/m^2)
262	!$OMP THREADPRIVATE(delsurfstor)
263
264	!
265	! Specific variables for simple routing
266	!
267	REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid)
268	!$OMP THREADPRIVATE(topoind_r)
269	INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo)
270	!$OMP THREADPRIVATE(route_flow_rp1)
271	REAL(r_std),SAVE,ALLOCATABLE :: fast_reservoir_r(:) !! Water amount in the fast reservoir (kg) - (local routing grid)
272	!$OMP THREADPRIVATE(fast_reservoir_r)
273	REAL(r_std),SAVE,ALLOCATABLE :: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid)
274	!$OMP THREADPRIVATE(slow_reservoir_r)
275	REAL(r_std),SAVE,ALLOCATABLE :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
276	!$OMP THREADPRIVATE(stream_reservoir_r)
277	LOGICAL,SAVE,ALLOCATABLE :: is_lakeinflow_r(:) !! is lake inflow point - (local routing grid)
278	!$OMP THREADPRIVATE(is_lakeinflow_r)
279	LOGICAL,SAVE,ALLOCATABLE :: is_coastalflow_r(:) !! is coastal flow point - (local routing grid)
280	!$OMP THREADPRIVATE(is_coastalflow_r)
281	LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid)
282	!$OMP THREADPRIVATE(is_riverflow_r)
283	LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid)
284	!$OMP THREADPRIVATE(is_streamflow_r)
285	LOGICAL,SAVE,ALLOCATABLE :: routing_mask_r(:) !! valid routing point - (local routing grid)
286	!$OMP THREADPRIVATE(routing_mask_r)
287	LOGICAL,SAVE,ALLOCATABLE :: coast_mask(:) !! is a coast point - (local native grid)
288	!$OMP THREADPRIVATE(coast_mask)
289	INTEGER,SAVE :: nbpt_r !! number of point in local routing grid
290	!$OMP THREADPRIVATE(nbpt_r)
291	INTEGER,SAVE :: nbpt_rp1 !! number of point in local routing grid with halo of 1
292	!$OMP THREADPRIVATE(nbpt_rp1)
293	REAL(r_std),SAVE,ALLOCATABLE :: routing_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
294	!$OMP THREADPRIVATE(routing_weight)
295	REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_r(:) !! Weight to transfer lost water into intersecting cell when doing interpolation from local routing grid to local native grid (lakeinflow)
296	!$OMP THREADPRIVATE(routing_weight_r)
297	REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_coast_r(:) !! Weight to transfer lost water into intersecting cell on coast
298	!$OMP THREADPRIVATE(routing_weight_coast_r)
299	!! when doing interpolation from local routing grid to local native grid (river flow+coastal flow)
300	REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:) !! basins riverflow id (local routing grid)
301	!$OMP THREADPRIVATE(basins_extended_r)
302	INTEGER(i_std) :: basins_count !! number of basins (local routing grid)
303	!$OMP THREADPRIVATE(basins_count)
304	INTEGER(i_std) :: basins_out !! number of basins to output for diag
305	!$OMP THREADPRIVATE(basins_out)
306
307	INTEGER(i_std) :: split_routing !! time spliting for routing
308	!$OMP THREADPRIVATE(split_routing)
309
310
311	INTEGER(i_std), PARAMETER :: nb_stations=14
312	REAL(r_std),PARAMETER :: station_lon(nb_stations) = &
313	(/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000, 28.7167, &
314	15.3000, 66.5300, 89.6700, 86.5000, 127.6500, 3.3833, 117.6200 /)
315	REAL(r_std),PARAMETER :: station_lat(nb_stations) = &
316	(/ 61.9340, 32.3150, -1.9470, -5.1281, 67.4583, 8.1500, 45.2167, &
317	-4.3000, 66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /)
318	CHARACTER(LEN=17),PARAMETER :: station_name(nb_stations) = &
319	(/ "Pilot station ", "Vicksburg ", "Obidos ", &
320	"Itupiranga ", "Arctic red river ", "Puente Angostura ", &
321	"Ceatal Izmail ", "Kinshasa ", "Salekhard ", &
322	"Bahadurabad ", "Igarka ", "Kusur ", &
323	"Malanville ", "Datong " /)
324
325	CONTAINS
326
327	!! =============================================================================================================================
328	!! SUBROUTINE: routing_simple_xios_initialize
329	!!
330	!>\BRIEF Initialize xios dependant defintion before closing context defintion
331	!!
332	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion.
333	!! This subroutine is called before the xios context is closed.
334	!!
335	!! RECENT CHANGE(S): None
336	!!
337	!! REFERENCE(S): None
338	!!
339	!! FLOWCHART: None
340	!! \n
341	!_ ==============================================================================================================================
342
343	SUBROUTINE routing_simple_xios_initialize
344	USE xios
345	USE routing, ONLY : routing_names
346	IMPLICIT NONE
347
348	!! 0 Variable and parameter description
349	CHARACTER(LEN=60),ALLOCATABLE :: label(:)
350	LOGICAL :: file_exists
351
352	IF (is_omp_root) THEN
353	CALL xios_get_axis_attr("basins", n_glo=basins_out) ! get nb basins to output
354	ALLOCATE(label(basins_out))
355	CALL routing_names(basins_out,label)
356	CALL xios_set_axis_attr("basins", label=label) ! set riverflow basins name
357	INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
358	IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.)
359	ENDIF
360
361	END SUBROUTINE routing_simple_xios_initialize
362
363
364
365	!! =============================================================================================================================
366	!! SUBROUTINE: routing_simple_init_2
367	!!
368	!>\BRIEF Privat subroutine
369	!!
370	!! DESCRIPTION: Privat subroutine to the module routing_simple. This subroutine is called in the end of routing_simple_initialize
371	!!
372	!! RECENT CHANGE(S): None
373	!!
374	!! REFERENCE(S): None
375	!!
376	!! FLOWCHART: None
377	!! \n
378	!_ ==============================================================================================================================
379
380	SUBROUTINE routing_simple_init_2(nbpt, contfrac)
381	USE xios
382	USE grid, ONLY : area
383	IMPLICIT NONE
384	INCLUDE "mpif.h"
385
386	!! 0 Variable and parameter description
387	!! 0.1 Input variables
388	INTEGER,INTENT(in) :: nbpt !! nb points native grid
389	REAL,INTENT(in) :: contfrac(nbpt) !! fraction of land
390
391	!! 0.2 Local variables
392	INTEGER :: ni !! longitude dimension of local routing grid
393	INTEGER :: nj !! latitude dimension of local routing grid
394
395	REAL(r_std),ALLOCATABLE :: trip_rp1(:,: ) !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid + halo (0:ni+1,0:nj+1)
396	REAL(r_std),ALLOCATABLE :: trip_extended_r(:,:) !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid (ni,nj)
397	!! routing is artificially computed on sea and endoric basins
398	LOGICAL,ALLOCATABLE :: routing_mask_rp1(:,:) !! valid routing point - local routing grid+halo (0:ni+1,0:nj+1)
399	REAL(r_std),ALLOCATABLE :: mask_native(:) !! some mask on native grid (nbpt)
400	REAL(r_std),ALLOCATABLE :: frac_routing_r(:) !! fraction of routing cell intesected by cells of native grid
401	REAL(r_std),ALLOCATABLE :: frac_routing_coast_r(:) !! fraction of routing cell intesected by coastal cells of native grid
402	INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir
403	REAL(r_std) :: sum1,sum2,sum_frac_routing_r
404	REAL(r_std) :: contfrac_mpi(nbp_mpi)
405	REAL(r_std) :: area_mpi(nbp_mpi)
406	INTEGER :: basins_count_mpi
407	INTEGER :: ierr
408	LOGICAL :: file_exists
409
410	!_ ================================================================================================================================
411
412	split_routing=1
413	CALL getin_p("SPLIT_ROUTING",split_routing)
414
415	CALL gather_omp(contfrac,contfrac_mpi)
416	CALL gather_omp(area, area_mpi)
417
418	IF (is_omp_root) THEN
419
420	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension
421
422	nbpt_r= ni*nj
423	nbpt_rp1= (ni+2)*(nj+2)
424
425	ALLOCATE(trip_extended_r(ni,nj))
426	ALLOCATE(topoind_r(ni*nj))
427	ALLOCATE(is_lakeinflow_r(ni*nj))
428	ALLOCATE(is_coastalflow_r(ni*nj))
429	ALLOCATE(is_riverflow_r(ni*nj))
430	ALLOCATE(is_streamflow_r(ni*nj))
431	ALLOCATE(routing_mask_r(ni*nj))
432	ALLOCATE(fast_reservoir_r(ni*nj))
433	ALLOCATE(slow_reservoir_r(ni*nj))
434	ALLOCATE(stream_reservoir_r(ni*nj))
435
436	ALLOCATE(mask_native(nbp_mpi))
437	ALLOCATE(coast_mask(nbp_mpi))
438	ALLOCATE(frac_routing_r(nbpt_r))
439	ALLOCATE(frac_routing_coast_r(nbpt_r))
440	ALLOCATE(routing_weight(nbp_mpi))
441	ALLOCATE(routing_weight_r(nbpt_r))
442	ALLOCATE(routing_weight_coast_r(nbpt_r))
443	ALLOCATE(basins_extended_r(nbpt_r))
444
445	coast_mask=.FALSE.
446	WHERE( contfrac_mpi(:)< 1.-1.e-5) coast_mask(:)=.TRUE. ! create mask for coastal cells on native grid
447
448	INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
449
450	IF (file_exists) THEN
451	CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r)
452	CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r)
453	CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r)
454	ELSE
455	fast_reservoir_r(:)=0
456	slow_reservoir_r(:)=0
457	stream_reservoir_r(:)=0
458	ENDIF
459
460	ALLOCATE(trip_rp1(0:ni+1,0:nj+1))
461	trip_rp1(:,:)=1e10
462	ALLOCATE(routing_mask_rp1(0:ni+1,0:nj+1))
463	ALLOCATE(route_flow_rp1((ni+2)*(nj+2)))
464
465	CALL xios_send_field("routing_contfrac",contfrac_mpi) ! diag
466	CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj)) ! recv trip array with halo of 1
467	CALL xios_recv_field("trip_extended_r",trip_extended_r) ! recv extended trip array from file
468	CALL xios_recv_field("topoind_r",topoind_r) ! recv topo index array from file
469	CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file
470
471	is_lakeinflow_r(:) = .FALSE.
472	is_coastalflow_r(:) = .FALSE.
473	is_riverflow_r(:) = .FALSE.
474
475	route_flow_rp1(:)=-1
476
477	routing_mask_rp1=.TRUE.
478	DO j=1,nj
479	DO i=1,ni
480	IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE. ! if ocean or endoric basins => masked point
481	ENDDO
482	ENDDO
483
484	mask_native(:)=0
485	frac_routing_r(:)=0
486	WHERE( .NOT. coast_mask(:)) mask_native(:)=1
487	CALL xios_send_field("mask_native_lake",mask_native) ! send full land point to XIOS (native grid)
488	CALL xios_recv_field("frac_routing_lake_r",frac_routing_r) ! receive fraction of intersected cell by full land, on routing grid
489
490
491	DO j=1,nj
492	DO i=1,ni
493	ij=i+ni*(j-1)
494	IF (frac_routing_r(ij)>1-1e-5) THEN ! if full land point (no part of coast)
495	IF ((.NOT. routing_mask_rp1(i,j)) .OR. trip_rp1(i,j)==98 .OR. trip_rp1(i,j)==99) THEN ! if masked point or river flow or coastal flow
496	routing_mask_rp1(i,j)=.TRUE. ! transform to non masked
497	trip_rp1(i,j)=trip_extended_r(i,j) ! replace by extended trip value
498	ENDIF
499	ENDIF
500	ENDDO
501	ENDDO
502
503	mask_native(:)=0
504	frac_routing_coast_r(:)=0
505	WHERE( coast_mask(:)) mask_native(:)=1
506	CALL xios_send_field("mask_native_coast",mask_native) ! send only coastal point to XIOS (native grid)
507	CALL xios_recv_field("frac_routing_coast_r",frac_routing_coast_r) ! receive fraction of intersected cell by coastal point, on routing grid
508
509	DO j=1,nj
510	DO i=1,ni
511	ij=i+ni*(j-1)
512	IF (frac_routing_coast_r(ij) > 0) THEN ! fraction of coastal point ?
513	IF (.NOT. routing_mask_rp1(i,j)) THEN ! if masked, transform to coastal flow to conserve the water balance
514	routing_mask_rp1(i,j)=.TRUE.
515	trip_rp1(i,j)=98 ! transform to coastal flow
516	topoind_r(ij)=1e-10 ! no residence time
517	ENDIF
518	ENDIF
519	ENDDO
520	ENDDO
521
522
523	mask_native(:)=1
524	frac_routing_r(:)=0
525	CALL xios_send_field("mask_native",mask_native) ! send only all point to XIOS (native grid)
526	CALL xios_recv_field("frac_routing_r",frac_routing_r) ! receive fraction of intersected cell, (routing grid)
527	DO j=1,nj
528	DO i=1,ni
529	ij=i+ni*(j-1)
530	IF (frac_routing_r(ij)>0) THEN
531	IF (.NOT. routing_mask_rp1(i,j)) THEN
532	routing_mask_rp1(i,j)=.TRUE. ! may never happen
533	ENDIF
534	ELSE
535	routing_mask_rp1(i,j)=.FALSE. ! if no intersection from native cells, mask the point (no routage)
536	trip_rp1(i,j)=1e10
537	ENDIF
538	ENDDO
539	ENDDO
540
541	CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo
542	CALL xios_recv_field("trip_rp1",trip_rp1) ! recv trip array with halo of 1
543
544
545	routing_mask_rp1=.TRUE.
546	DO j=0,nj+1
547	DO i=0,ni+1
548	IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE. ! update mask with halo of 1
549	ENDDO
550	ENDDO
551
552	!! Compute the routing
553	!! Loop on all point point of the local routing grid + halo
554	DO j=0,nj+1
555	jp1=j+1
556	jm1=j-1
557	DO i=0,ni+1
558	ij_rp1=i+(ni+2)*j+1
559	ij_r=i+ni*(j-1)
560	ip1=i+1
561	im1=i-1
562
563	IF (trip_rp1(i,j) < 100) THEN
564
565	ir=-1 ; jr=-1 ! -> -1 for 97,98,99
566	SELECT CASE (NINT(trip_rp1(i,j))) ! get the trip value for each points
567	CASE (1)
568	jr=jm1 ; ir=i ! north
569	CASE (2)
570	jr=jm1 ; ir=ip1 ! north-east
571	CASE (3)
572	jr=j ; ir=ip1 ! east
573	CASE (4)
574	jr=jp1 ; ir=ip1 ! south-east
575	CASE (5)
576	jr=jp1 ; ir=i ! south
577	CASE(6)
578	jr=jp1 ; ir=im1 ! south-west
579	CASE (7)
580	jr=j ; ir=im1 ! west
581	CASE (8)
582	jr=jm1 ; ir=im1 ! north-west
583	CASE (97)
584	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
585	is_lakeinflow_r(ij_r)=.TRUE. ! I am a lakeinflow point and route to myself
586	jr=j ; ir=i
587	ENDIF
588	CASE (98)
589	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
590	is_coastalflow_r(ij_r)=.TRUE. ! I am a coastal flow point and route to myself
591	jr=j ; ir=i
592	ENDIF
593	CASE (99)
594	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
595	is_riverflow_r(ij_r)=.TRUE. ! I am a riverflow point and route to myself
596	jr=j ; ir=i
597	ENDIF
598	END SELECT
599
600	IF (ir<0 .OR. ir>ni+1 .OR. jr<0 .OR. jr>nj+1) THEN
601	route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain+halo, no routing (will be done by other process)
602	ELSE
603
604	IF ( .NOT. routing_mask_rp1(ir,jr)) THEN ! if routing to a masked point, cut the flow
605	jr=j ; ir=i ; ! and route to myself
606	IF (i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if some flow come from neighbour cell
607	IF ( trip_rp1(i-1,j-1)==4 .OR. trip_rp1(i,j-1)==5 .OR. trip_rp1(i+1,j)==6 &
608	.OR. trip_rp1(i-1,j)==3 .OR. trip_rp1(i+1,j)==7 &
609	.OR. trip_rp1(i-1,j+1)==2 .OR. trip_rp1(i-1,j+1)==1 .OR. trip_rp1(i+1,j+1)==8) THEN
610	is_riverflow_r(ij_r)=.TRUE. ! => transform to riverflow
611	ELSE
612	is_coastalflow_r(ij_r)=.TRUE. ! else transform to coastalflow
613	ENDIF
614	ENDIF
615	ENDIF
616
617	IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN
618	route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain, no routing (will be done by other process)
619	ELSE
620	route_flow_rp1(ij_rp1)=ir+ni*(jr-1) ! define the cell where to flow
621	ENDIF
622	ENDIF
623
624	ENDIF
625	ENDDO
626	ENDDO
627
628	routing_mask_r(:)=reshape(routing_mask_rp1(1:ni,1:nj),(/ni*nj/))
629
630	is_streamflow_r(:)=.NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:)
631
632	DO ij=1,nbp_mpi
633	routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij)
634	ENDDO
635
636	routing_weight_r(:)=0
637	DO ij=1,nbpt_r
638	IF (frac_routing_r(ij)>0) routing_weight_r(ij)=1./frac_routing_r(ij)
639	ENDDO
640
641	routing_weight_coast_r(:)=0.
642	DO ij=1,nbpt_r
643	IF (frac_routing_coast_r(ij)>0) routing_weight_coast_r(ij)=1./frac_routing_coast_r(ij)
644	ENDDO
645
646	CALL xios_send_field("routing_weight_coast_r",routing_weight_coast_r)
647	! looking for basins
648	basins_count_mpi=0
649	DO ij=1,nbpt_r
650	IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij)
651	ENDDO
652	CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr)
653
654	ELSE
655
656	nbpt_r=0
657
658	ENDIF
659
660	END SUBROUTINE routing_simple_init_2
661
662
663	!! ================================================================================================================================
664	!! SUBROUTINE : routing_simple_flow
665	!!
666	!>\BRIEF This subroutine computes the transport of water in the various reservoirs
667	!! (including ponds and floodplains) and the water withdrawals from the reservoirs for irrigation.
668	!!
669	!! DESCRIPTION (definitions, functional, design, flags) :
670	!! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an
671	!! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes
672	!! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of
673	!! the up-stream basin and adding it to the down-stream one.
674	!! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once
675	!! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them
676	!! the hydrographs of the largest rivers we have chosen to monitor.
677	!!
678	!! RECENT CHANGE(S): None
679	!!
680	!! MAIN OUTPUT VARIABLE(S): lakeinflow, returnflow, reinfiltration, irrigation, riverflow, coastalflow, hydrographs, flood_frac, flood_res
681	!!
682	!! REFERENCES :
683	!! - Ngo-Duc, T., K. Laval, G. Ramillien, J. Polcher, and A. Cazenave (2007)
684	!! Validation of the land water storage simulated by Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) with Gravity Recovery and Climate Experiment (GRACE) data.
685	!! Water Resour. Res., 43, W04427, doi:10.1029/2006WR004941.
686	!! * Irrigation:
687	!! - de Rosnay, P., J. Polcher, K. Laval, and M. Sabre (2003)
688	!! Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula.
689	!! Geophys. Res. Lett., 30(19), 1986, doi:10.1029/2003GL018024.
690	!! - A.C. Vivant (2003)
691	!! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte.
692	!! , , 51pp.
693	!! - N. Culson (2004)
694	!! Impact de l'irrigation sur le cycle de l'eau
695	!! Master thesis, Paris VI University, 55pp.
696	!! - X.-T. Nguyen-Vinh (2005)
697	!! Analyse de l'impact de l'irrigation en Amerique du Nord - plaine du Mississippi - sur la climatologie regionale
698	!! Master thesis, Paris VI University, 33pp.
699	!! - M. Guimberteau (2006)
700	!! Analyse et modifications proposees de la modelisation de l'irrigation dans un modele de surface.
701	!! Master thesis, Paris VI University, 46pp.
702	!! - Guimberteau M. (2010)
703	!! Modelisation de l'hydrologie continentale et influences de l'irrigation sur le cycle de l'eau.
704	!! Ph.D. thesis, Paris VI University, 195pp.
705	!! - Guimberteau M., Laval K., Perrier A. and Polcher J. (2011).
706	!! Global effect of irrigation and its impact on the onset of the Indian summer monsoon.
707	!! In press, Climate Dynamics, doi: 10.1007/s00382-011-1252-5.
708	!! * Floodplains:
709	!! - A.C. Vivant (2002)
710	!! L'ecoulement lateral de l'eau sur les surfaces continentales. Prise en compte des plaines d'inondations dans ORCHIDEE.
711	!! Master thesis, Paris VI University, 46pp.
712	!! - A.C. Vivant (2003)
713	!! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte.
714	!! , , 51pp.
715	!! - T. d'Orgeval (2006)
716	!! Impact du changement climatique sur le cycle de l'eau en Afrique de l'Ouest: modelisation et incertitudes.
717	!! Ph.D. thesis, Paris VI University, 188pp.
718	!! - T. d'Orgeval, J. Polcher, and P. de Rosnay (2008)
719	!! Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes.
720	!! Hydrol. Earth Syst. Sci., 12, 1387-1401
721	!! - M. Guimberteau, G. Drapeau, J. Ronchail, B. Sultan, J. Polcher, J.-M. Martinez, C. Prigent, J.-L. Guyot, G. Cochonneau,
722	!! J. C. Espinoza, N. Filizola, P. Fraizy, W. Lavado, E. De Oliveira, R. Pombosa, L. Noriega, and P. Vauchel (2011)
723	!! Discharge simulation in the sub-basins of the Amazon using ORCHIDEE forced by new datasets.
724	!! Hydrol. Earth Syst. Sci. Discuss., 8, 11171-11232, doi:10.5194/hessd-8-11171-2011
725	!!
726	!! FLOWCHART :None
727	!! \n
728	!_ ================================================================================================================================
729
730	SUBROUTINE routing_simple_flow(nbpt, dt_routing, lalo, floodout, runoff_omp, drainage_omp, &
731	vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, &
732	lakeinflow_omp, returnflow, reinfiltration, irrigation, riverflow_omp, &
733	coastalflow_omp, hydrographs, slowflow_diag, flood_frac, flood_res)
734
735	USE xios
736	USE grid, ONLY : area
737	IMPLICIT NONE
738	INCLUDE "mpif.h"
739
740	!! 0 Variable and parameter description
741	!! 0.1 Input variables
742
743	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
744	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
745	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes
746	REAL(r_std), INTENT(in) :: runoff_omp(nbpt) !! Grid-point runoff (kg/m^2/dt)
747	REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt)
748	REAL(r_std), INTENT(in) :: drainage_omp(nbpt) !! Grid-point drainage (kg/m^2/dt)
749	REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated fraction (unitless;0-1)
750	REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas which can not have vegetation
751	REAL(r_std), INTENT(in) :: transpot_mean(nbpt) !! Mean potential transpiration of the vegetation (kg/m^2/dt)
752	REAL(r_std), INTENT(in) :: precip(nbpt) !! Rainfall (kg/m^2/dt)
753	REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless)
754	REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
755	REAL(r_std), INTENT(in) :: floodtemp(nbpt) !! Temperature to decide if floodplains work (K)
756	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)
757
758	!! 0.2 Output variables
759
760	REAL(r_std), INTENT(out) :: lakeinflow_omp(nbpt) !! Water inflow to the lakes (kg/dt)
761	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns into the grid box.
762	!! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
763	REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
764	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)
765	REAL(r_std), INTENT(out) :: riverflow_omp(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
766	REAL(r_std), INTENT(out) :: coastalflow_omp(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
767	REAL(r_std), INTENT(out) :: hydrographs(nbpt) !! Hydrographs at the outflow of the grid box for major basins (kg/dt)
768	REAL(r_std), INTENT(out) :: slowflow_diag(nbpt) !! Hydrographs of slow_flow = routed slow_flow for major basins (kg/dt)
769	REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
770	REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
771
772	!! 0.4 Local variables
773	REAL(r_std) :: runoff(nbp_mpi) !! Grid-point runoff (kg/dt)
774	REAL(r_std) :: drainage(nbp_mpi) !! Grid-point drainage (kg/dt)
775	REAL(r_std) :: riverflow(nbp_mpi)
776	REAL(r_std) :: coastalflow(nbp_mpi)
777	REAL(r_std) :: lakeinflow(nbp_mpi)
778	REAL(r_std) :: fast_diag_mpi(nbp_mpi)
779	REAL(r_std) :: slow_diag_mpi(nbp_mpi)
780	REAL(r_std) :: stream_diag_mpi(nbp_mpi)
781
782	! from input model -> routing_grid
783	REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt)
784	REAL(r_std) :: floodout_r(nbpt_r) !! Grid-point flow out of floodplains (kg/m^2/dt)
785	REAL(r_std) :: drainage_r(nbpt_r) !! Grid-point drainage (kg/m^2/dt)
786	REAL(r_std) :: vegtot_r(nbpt_r) !! Potentially vegetated fraction (unitless;0-1)
787	REAL(r_std) :: totnobio_r(nbpt_r) !! Other areas which can not have vegetation
788	REAL(r_std) :: transpot_mean_r(nbpt_r) !! Mean potential transpiration of the vegetation (kg/m^2/dt)
789	REAL(r_std) :: precip_r(nbpt_r) !! Rainfall (kg/m^2/dt)
790	REAL(r_std) :: humrel_r(nbpt_r) !! Soil moisture stress, root extraction potential (unitless)
791	REAL(r_std) :: k_litt_r(nbpt_r) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
792	REAL(r_std) :: floodtemp_r(nbpt_r) !! Temperature to decide if floodplains work (K)
793	REAL(r_std) :: reinf_slope_r(nbpt_r) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
794
795
796
797	REAL(r_std), DIMENSION(nbpt_r) :: fast_flow_r !! Outflow from the fast reservoir (kg/dt)
798	REAL(r_std), DIMENSION(nbpt_r) :: slow_flow_r !! Outflow from the slow reservoir (kg/dt)
799	REAL(r_std), DIMENSION(nbpt_r) :: stream_flow_r !! Outflow from the stream reservoir (kg/dt)
800	REAL(r_std), DIMENSION(nbpt_r) :: hydrographs_r !! hydrograph (kg/dt)
801	REAL(r_std), DIMENSION(nbpt_r) :: flood_flow_r !! Outflow from the floodplain reservoir (kg/dt)
802	REAL(r_std), DIMENSION(nbpt_r) :: pond_inflow_r !! Inflow to the pond reservoir (kg/dt)
803	REAL(r_std), DIMENSION(nbpt_r) :: pond_drainage_r !! Drainage from pond (kg/m^2/dt)
804	REAL(r_std), DIMENSION(nbpt_r) :: flood_drainage_r !! Drainage from floodplains (kg/m^2/dt)
805	REAL(r_std), DIMENSION(nbpt_r) :: return_swamp_r !! Inflow to the swamp (kg/dt)
806	!
807	! Irrigation per basin
808	!
809	REAL(r_std), DIMENSION(nbpt_r) :: irrig_needs_r !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
810	REAL(r_std), DIMENSION(nbpt_r) :: irrig_actual_r !! Possible irrigation according to the water availability in the reservoirs (kg)
811	REAL(r_std), DIMENSION(nbpt_r) :: irrig_deficit_r !! Amount of water missing for irrigation (kg)
812	REAL(r_std), DIMENSION(nbpt_r) :: irrig_adduct_r !! Amount of water carried over from other basins for irrigation (kg)
813	!
814	REAL(r_std), DIMENSION(nbpt_r) :: transport_r !! Water transport between basins (kg/dt)
815	REAL(r_std), DIMENSION(nbpt_r) :: floods_r !! Water flow in to the floodplains (kg/dt)
816	REAL(r_std), DIMENSION(nbpt_r) :: potflood_r !! Potential inflow to the swamps (kg/dt)
817	REAL(r_std), DIMENSION(nbpt) :: tobeflooded !! Maximal surface which can be inundated in each grid box (m^2)
818	REAL(r_std), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2)
819	REAL(r_std), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg) (not used)
820	REAL(r_std), DIMENSION(nbasmax) :: pond_excessflow !!
821	REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt)
822	REAL(r_std) :: floodindex !! Fraction of grid box area inundated (unitless;0-1)
823	REAL(r_std) :: pondex !!
824	REAL(r_std) :: flood_frac_pot !! Total fraction of the grid box which is flooded at optimum repartition (unitless;0-1)
825	REAL(r_std) :: stream_tot !! Total water amount in the stream reservoirs (kg)
826	REAL(r_std) :: adduction !! Importation of water from a stream reservoir of a neighboring grid box (kg)
827	REAL(r_std), DIMENSION(8,nbasmax) :: streams_around !! Stream reservoirs of the neighboring grid boxes (kg)
828	INTEGER(i_std), DIMENSION(8) :: igrd !!
829	INTEGER(i_std), DIMENSION(2) :: ff !!
830	INTEGER(i_std), DIMENSION(1) :: fi !!
831	INTEGER(i_std) :: ig, ib, ib2, ig2 !! Indices (unitless)
832	INTEGER(i_std) :: rtg, rtb, in !! Indices (unitless)
833	INTEGER(i_std) :: ier !! Error handling
834	INTEGER(i_std) :: isplit
835
836	LOGICAL, PARAMETER :: check_reservoir = .TRUE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false)
837	REAL(r_std), DIMENSION(nbpt_r) :: flood_frac_bas_r
838	REAL(r_std), DIMENSION(nbpt_r) :: flood_reservoir_r
839	REAL(r_std), DIMENSION(nbpt_r) :: pond_reservoir_r
840
841	REAL(r_std), DIMENSION(nbpt_r) :: lakeinflow_r
842	REAL(r_std), DIMENSION(nbpt_r) :: coastalflow_r
843	REAL(r_std), DIMENSION(nbpt_r) :: riverflow_r
844
845	REAL(r_std), DIMENSION(nbpt_r) :: flow_r
846	REAL(r_std), DIMENSION(nbpt_rp1) :: flow_rp1
847	REAL(r_std) :: water_balance_before, water_balance_after
848	REAL(r_std) :: basins_riverflow_mpi(0:basins_count)
849	REAL(r_std) :: basins_riverflow(0:basins_count)
850
851	INTEGER :: ierr
852
853	!_ ================================================================================================================================
854	!
855
856
857	irrig_netereq(:) = zero
858	totarea(:) = zero
859	totflood(:) = zero
860
861	flood_frac_bas_r(:)=zero
862	!
863	!> The outflow fluxes from the three reservoirs are computed.
864	!> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume.
865	!> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid
866	!> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula
867	!> (Dingman, 1994; Ducharne et al., 2003).
868	!> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day)
869	!> which is an e-folding time, the time necessary for the water amount
870	!> in the stream reservoir to decrease by a factor e. Hence, it gives an order of
871	!> magnitude of the travel time through this reservoir between
872	!> the sub-basin considered and its downstream neighbor.
873	CALL gather_omp(runoff_omp,runoff)
874	CALL gather_omp(drainage_omp, drainage)
875
876	IF (is_omp_root) THEN
877
878	hydrographs_r(:)=0
879
880	DO isplit=1,split_routing
881
882
883	irrig_needs_r(:) = zero
884	irrig_actual_r(:) = zero
885	irrig_deficit_r(:) = zero
886	irrig_adduct_r(:) = zero
887
888	DO ig=1,nbpt_r
889	IF ( routing_mask_r(ig) ) THEN
890	!
891	! Each of the fluxes is limited by the water in the reservoir and a small margin
892	! (min_reservoir) to avoid rounding errors.
893	!
894	! IF (fast_reservoir_r(ig)>1e-30) THEN
895	! PRINT*,"fast_reservoir > 0 ",ig, fast_reservoir_r(ig), topoind_r(ig)
896	! ENDIF
897	flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)fast_tcstone_day/(dt_routing/split_routing)),&
898	& fast_reservoir_r(ig)-min_sechiba)
899	fast_flow_r(ig) = MAX(flow, zero)
900	!
901	flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)slow_tcstone_day/(dt_routing/split_routing)),&
902	& slow_reservoir_r(ig)-min_sechiba)
903	slow_flow_r(ig) = MAX(flow, zero)
904	!
905	flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)stream_tcst &
906	& MAX(un-SQRT(flood_frac_bas_r(ig)),min_sechiba)*one_day/(dt_routing/split_routing)),&
907	& stream_reservoir_r(ig)-min_sechiba)
908	stream_flow_r(ig) = MAX(flow, zero)
909	!
910	!
911	ELSE
912	fast_flow_r(ig) = zero
913	slow_flow_r(ig) = zero
914	stream_flow_r(ig) = zero
915	ENDIF
916	ENDDO
917
918	DO ig=1,nbpt_r
919	flood_drainage_r(ig) = zero
920	flood_flow_r(ig) = zero
921	flood_reservoir_r(ig) = zero
922	ENDDO
923
924	DO ig=1,nbpt_r
925	pond_inflow_r(ig) = zero
926	pond_drainage_r(ig) = zero
927	pond_reservoir_r(ig) = zero
928	ENDDO
929
930
931	!-
932	!- Compute the transport
933	!-
934
935	flow_r(:)=fast_flow_r(:) + slow_flow_r(:) + stream_flow_r(:)
936
937	CALL xios_send_field("flow_r",flow_r) ! transfer halo
938	CALL xios_recv_field("flow_rp1",flow_rp1)
939
940	transport_r(:)=0
941
942	DO ig=1,nbpt_rp1
943	IF ( route_flow_rp1(ig) > 0 ) THEN
944	transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig)
945	ENDIF
946	ENDDO
947
948
949	!-
950	!- Do the floodings - First initialize
951	!-
952	return_swamp_r(:)=zero
953	floods_r(:)=zero
954
955	!
956	! Update all reservoirs
957	!> The slow and deep reservoir (slow_reservoir) collect the deep drainage whereas the
958	!> fast_reservoir collects the computed surface runoff. Both discharge into a third reservoir
959	!> (stream_reservoir) of the next sub-basin downstream.
960	!> Water from the floodplains reservoir (flood_reservoir) flows also into the stream_reservoir of the next sub-basin downstream.
961	!> Water that flows into the pond_reservoir is withdrawn from the fast_reservoir.
962	!
963	runoff=runoff*routing_weight
964	CALL xios_send_field("routing_runoff",runoff) ! interp conservative model -> routing
965	CALL xios_recv_field("routing_runoff_r",runoff_r)
966	drainage=drainage*routing_weight
967	CALL xios_send_field("routing_drainage",drainage) ! interp conservative model -> routing
968	CALL xios_recv_field("routing_drainage_r",drainage_r)
969
970	CALL xios_send_field("fast_reservoir_r",fast_reservoir_r)
971	CALL xios_send_field("slow_reservoir_r",slow_reservoir_r)
972	CALL xios_send_field("stream_reservoir_r",stream_reservoir_r)
973
974	WHERE (.NOT. routing_mask_r) runoff_r(:)=0
975	WHERE (.NOT. routing_mask_r) drainage_r(:)=0
976
977	CALL MPI_ALLREDUCE(sum(runoff+drainage),water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
978	CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
979
980	PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, &
981	" ; delta : ", 100.(water_balance_after-water_balance_before)/(0.5(water_balance_after+water_balance_before)),"%"
982
983	CALL xios_recv_field("water_balance_before",water_balance_before)
984
985
986	DO ig=1,nbpt_r
987	IF ( routing_mask_r(ig) ) THEN
988
989	!
990	fast_reservoir_r(ig) = fast_reservoir_r(ig) + runoff_r(ig) - &
991	& fast_flow_r(ig) - pond_inflow_r(ig)
992	!
993	slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - &
994	& slow_flow_r(ig)
995	!
996
997	stream_reservoir_r(ig) = stream_reservoir_r(ig) + flood_flow_r(ig) - &
998	& stream_flow_r(ig) - return_swamp_r(ig) - floods_r(ig)
999
1000	IF (is_streamflow_r(ig)) stream_reservoir_r(ig)= stream_reservoir_r(ig) + transport_r(ig)
1001	!
1002	flood_reservoir_r(ig) = flood_reservoir_r(ig) + floods_r(ig) - &
1003	& flood_flow_r(ig)
1004	!
1005	pond_reservoir_r(ig) = pond_reservoir_r(ig) + pond_inflow_r(ig) - pond_drainage_r(ig)
1006	!
1007	IF ( flood_reservoir_r(ig) .LT. zero ) THEN
1008	IF ( check_reservoir ) THEN
1009	WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ". Problem is being corrected."
1010	WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir_r(ig), floods_r(ig), &
1011	& flood_flow_r(ig)
1012	ENDIF
1013	stream_reservoir_r(ig) = stream_reservoir_r(ib) + flood_reservoir_r(ib)
1014	flood_reservoir_r(ib) = zero
1015	ENDIF
1016	!
1017	IF ( stream_reservoir_r(ig) .LT. zero ) THEN
1018	IF ( check_reservoir ) THEN
1019	WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected."
1020	WRITE(numout,*) "stream_reservoir, flood_flow, transport : ", stream_reservoir_r(ig), flood_flow_r(ig), &
1021	& transport_r(ig)
1022	WRITE(numout,*) "stream_flow, return_swamp, floods :", stream_flow_r(ig), return_swamp_r(ig), floods_r(ig)
1023	ENDIF
1024	fast_reservoir_r(ig) = fast_reservoir_r(ig) + stream_reservoir_r(ig)
1025	stream_reservoir_r(ig) = zero
1026	ENDIF
1027	!
1028	IF ( fast_reservoir_r(ig) .LT. zero ) THEN
1029	IF ( check_reservoir ) THEN
1030	WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected."
1031	WRITE(numout,*) "fast_reservoir, runoff, fast_flow, ponf_inflow : ", fast_reservoir_r(ig), &
1032	&runoff_r(ig), fast_flow_r(ig), pond_inflow_r(ig)
1033	ENDIF
1034	slow_reservoir_r(ig) = slow_reservoir_r(ig) + fast_reservoir_r(ig)
1035	fast_reservoir_r(ig) = zero
1036	ENDIF
1037
1038	IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN
1039	WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig
1040	WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', &
1041	& slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig)
1042	WRITE(numout,*) 'WARNING : pondr, pond_inflow, pond_drainage', &
1043	& pond_reservoir_r(ig), pond_inflow_r(ig), pond_drainage_r(ig)
1044	CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','')
1045	ENDIF
1046
1047	ENDIF
1048	ENDDO
1049
1050	DO ig=1,nbpt_r
1051	IF ( routing_mask_r(ig) ) THEN
1052	hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig)
1053	ENDIF
1054	ENDDO
1055
1056	ENDDO ! isplit
1057
1058	!
1059	!
1060	!
1061	! The three output types for the routing : endoheric basins,, rivers and
1062	! diffuse coastal flow.
1063	!
1064	lakeinflow_r(:)=0
1065	coastalflow_r(:)=0
1066	riverflow_r(:)=0
1067	basins_riverflow_mpi(:)=0
1068
1069	DO ig=1,nbpt_r
1070	IF ( routing_mask_r(ig) ) THEN
1071
1072	IF (is_lakeinflow_r(ig)) THEN
1073	lakeinflow_r(ig) = transport_r(ig)
1074	basins_riverflow_mpi(basins_extended_r(ig)) = &
1075	basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig)
1076	ENDIF
1077
1078	IF (is_coastalflow_r(ig)) THEN
1079	coastalflow_r(ig) = transport_r(ig)
1080	basins_riverflow_mpi(basins_extended_r(ig)) = &
1081	basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig)
1082	ENDIF
1083
1084	IF (is_riverflow_r(ig)) THEN
1085	riverflow_r(ig) = transport_r(ig)
1086	basins_riverflow_mpi(basins_extended_r(ig)) = &
1087	basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig)
1088	ENDIF
1089
1090	ENDIF
1091	ENDDO
1092
1093
1094	CALL MPI_ALLREDUCE(basins_riverflow_mpi,basins_riverflow,basins_count+1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1095	CALL xios_send_field("basins_riverflow",basins_riverflow(1:basins_out)/1000./dt_routing)
1096
1097	!
1098	flood_res = flood_diag + pond_diag
1099	!
1100	CALL xios_send_field("routing_lakeinflow_r" ,lakeinflow_r*routing_weight_r)
1101	CALL xios_recv_field("routing_lakeinflow" ,lakeinflow)
1102	CALL xios_send_field("routing_coastalflow_r" ,coastalflow_r*routing_weight_coast_r)
1103	CALL xios_recv_field("routing_coastalflow_temp" ,coastalflow)
1104	WHERE(.NOT. coast_mask(:)) coastalflow(:)=0
1105	CALL xios_send_field("routing_coastalflow" ,coastalflow)
1106
1107	CALL xios_send_field("routing_riverflow_r" ,riverflow_r*routing_weight_coast_r)
1108	CALL xios_recv_field("routing_riverflow_temp" ,riverflow)
1109	WHERE(.NOT. coast_mask(:)) riverflow(:)=0
1110	CALL xios_send_field("routing_riverflow" ,riverflow)
1111
1112	CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow)
1113	CALL xios_send_field("routing_hydrographs_r", hydrographs_r/1000./dt_routing)
1114
1115	! diag
1116	CALL xios_send_field("routing_fast_reservoir_r" , fast_reservoir_r)
1117	CALL xios_recv_field("routing_fast_reservoir" , fast_diag_mpi)
1118
1119	CALL xios_send_field("routing_slow_reservoir_r" , slow_reservoir_r)
1120	CALL xios_recv_field("routing_slow_reservoir" , slow_diag_mpi)
1121
1122	CALL xios_send_field("routing_stream_reservoir_r" , stream_reservoir_r)
1123	CALL xios_recv_field("routing_stream_reservoir" , stream_diag_mpi)
1124
1125	CALL xios_recv_field("water_balance_after" ,water_balance_after)
1126
1127	PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, &
1128	" ; delta : ", 100(water_balance_after-water_balance_before)/(0.5(water_balance_after+water_balance_before)),"%"
1129
1130
1131
1132	ENDIF ! is_omp_root
1133
1134	CALL scatter_omp(riverflow,riverflow_omp)
1135	CALL scatter_omp(coastalflow,coastalflow_omp)
1136	CALL scatter_omp(lakeinflow,lakeinflow_omp)
1137	CALL scatter_omp(fast_diag_mpi,fast_diag)
1138	CALL scatter_omp(slow_diag_mpi,slow_diag)
1139	CALL scatter_omp(stream_diag_mpi,stream_diag)
1140
1141	DO ig=1,nbpt
1142	fast_diag(ig)=fast_diag(ig)/area(ig)
1143	slow_diag(ig)=slow_diag(ig)/area(ig)
1144	stream_diag(ig)=stream_diag(ig)/area(ig)
1145	ENDDO
1146
1147
1148	flood_frac(:) = zero
1149	flood_height(:) = zero
1150	flood_frac_bas_r(:) = zero
1151
1152	returnflow(:) = zero
1153	reinfiltration(:) = zero
1154
1155	!
1156	!
1157	! Compute the fluxes which leave the routing scheme
1158	!
1159	! Lakeinflow is in Kg/dt
1160	! returnflow is in Kg/m^2/dt
1161	!
1162	delsurfstor(:) = -flood_diag(:)-pond_diag(:)-lake_diag(:)
1163	hydrographs(:) = zero
1164	slowflow_diag(:) = zero
1165	flood_diag(:) = zero
1166	pond_diag(:) = zero
1167	irrigation(:) = zero
1168
1169
1170	END SUBROUTINE routing_simple_flow
1171
1172
1173	!! =============================================================================================================================
1174	!! SUBROUTINE: routing_simple_initialize
1175	!!
1176	!>\BRIEF Initialize the routing_simple module
1177	!!
1178	!! DESCRIPTION: Initialize the routing_simple module. Read from restart file or read the routing.nc file to initialize the
1179	!! routing scheme.
1180	!!
1181	!! RECENT CHANGE(S)
1182	!!
1183	!! REFERENCE(S)
1184	!!
1185	!! FLOWCHART
1186	!! \n
1187	!_ ==============================================================================================================================
1188
1189	SUBROUTINE routing_simple_initialize( kjit, nbpt, index, &
1190	rest_id, hist_id, hist2_id, lalo, &
1191	neighbours, resolution, contfrac, stempdiag, &
1192	returnflow, reinfiltration, irrigation, riverflow, &
1193	coastalflow, flood_frac, flood_res )
1194
1195	IMPLICIT NONE
1196
1197	!! 0 Variable and parameter description
1198	!! 0.1 Input variables
1199	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1200	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
1201	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless)
1202	INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless)
1203	INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless)
1204	INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless)
1205	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
1206
1207	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point
1208	!! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
1209	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m)
1210	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
1211	REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
1212
1213	!! 0.2 Output variables
1214	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box.
1215	!! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
1216	REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
1217	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)
1218	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)
1219
1220	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)
1221	REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
1222	REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
1223
1224	!! 0.3 Local variables
1225	LOGICAL :: init_irrig !! Logical to initialize the irrigation (true/false)
1226	LOGICAL :: init_flood !! Logical to initialize the floodplains (true/false)
1227	LOGICAL :: init_swamp !! Logical to initialize the swamps (true/false)
1228
1229	!_ ================================================================================================================================
1230
1231	!
1232	! do initialisation
1233	!
1234	nbvmax = 440
1235	! Here we will allocate the memory and get the fixed fields from the restart file.
1236	! If the info is not found then we will compute the routing map.
1237	!
1238	CALL routing_simple_init_1 (kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
1239	riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id)
1240
1241	routing_area => routing_area_loc
1242	topo_resid => topo_resid_loc
1243	route_togrid => route_togrid_loc
1244	route_tobasin => route_tobasin_loc
1245	global_basinid => global_basinid_loc
1246	hydrodiag => hydrodiag_loc
1247
1248	! This routine computes the routing map if the route_togrid_glo is undefined. This means that the
1249	! map has not been initialized during the restart process..
1250	!
1251	!! Reads in the map of the basins and flow directions to construct the catchments of each grid box
1252	!
1253	IF ( COUNT(route_togrid_glo .GE. undef_int) .GT. 0 ) THEN
1254	! CALL routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac)
1255	ENDIF
1256	!
1257	! Do we have what we need if we want to do irrigation
1258	!! Initialisation of flags for irrigated land, flood plains and swamps
1259	!
1260	init_irrig = .FALSE.
1261	IF ( do_irrigation ) THEN
1262	IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE.
1263	END IF
1264
1265	init_flood = .FALSE.
1266	IF ( do_floodplains ) THEN
1267	IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) init_flood = .TRUE.
1268	END IF
1269
1270	init_swamp = .FALSE.
1271	IF ( doswamps ) THEN
1272	IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE.
1273	END IF
1274
1275	!! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree
1276	!! base data set to the resolution of the model.
1277
1278	IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN
1279	! CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, &
1280	! contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)
1281	ENDIF
1282
1283	IF ( do_irrigation ) THEN
1284	CALL xios_orchidee_send_field("irrigmap",irrigated)
1285
1286	WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated)
1287	IF ( .NOT. almaoutput ) THEN
1288	CALL histwrite_p(hist_id, 'irrigmap', 1, irrigated, nbpt, index)
1289	ELSE
1290	CALL histwrite_p(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index)
1291	ENDIF
1292	IF ( hist2_id > 0 ) THEN
1293	IF ( .NOT. almaoutput ) THEN
1294	CALL histwrite_p(hist2_id, 'irrigmap', 1, irrigated, nbpt, index)
1295	ELSE
1296	CALL histwrite_p(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index)
1297	ENDIF
1298	ENDIF
1299	ENDIF
1300
1301	IF ( do_floodplains ) THEN
1302	CALL xios_orchidee_send_field("floodmap",floodplains)
1303
1304	WRITE(numout,*) 'Verification : range of floodplains : ', MINVAL(floodplains), MAXVAL(floodplains)
1305	IF ( .NOT. almaoutput ) THEN
1306	CALL histwrite_p(hist_id, 'floodmap', 1, floodplains, nbpt, index)
1307	ELSE
1308	CALL histwrite_p(hist_id, 'FloodplainsMap', 1, floodplains, nbpt, index)
1309	ENDIF
1310	IF ( hist2_id > 0 ) THEN
1311	IF ( .NOT. almaoutput ) THEN
1312	CALL histwrite_p(hist2_id, 'floodmap', 1, floodplains, nbpt, index)
1313	ELSE
1314	CALL histwrite_p(hist2_id, 'FloodplainsMap', 1, floodplains, nbpt, index)
1315	ENDIF
1316	ENDIF
1317	ENDIF
1318
1319	IF ( doswamps ) THEN
1320	CALL xios_orchidee_send_field("swampmap",swamp)
1321
1322	WRITE(numout,*) 'Verification : range of swamp : ', MINVAL(swamp), MAXVAL(swamp)
1323	IF ( .NOT. almaoutput ) THEN
1324	CALL histwrite_p(hist_id, 'swampmap', 1, swamp, nbpt, index)
1325	ELSE
1326	CALL histwrite_p(hist_id, 'SwampMap', 1, swamp, nbpt, index)
1327	ENDIF
1328	IF ( hist2_id > 0 ) THEN
1329	IF ( .NOT. almaoutput ) THEN
1330	CALL histwrite_p(hist2_id, 'swampmap', 1, swamp, nbpt, index)
1331	ELSE
1332	CALL histwrite_p(hist2_id, 'SwampMap', 1, swamp, nbpt, index)
1333	ENDIF
1334	ENDIF
1335	ENDIF
1336
1337	!! This routine gives a diagnostic of the basins used.
1338	! CALL routing_diagnostic_p(nbpt, index, lalo, resolution, contfrac, hist_id, hist2_id)
1339
1340	CALL routing_simple_init_2(nbpt, contfrac)
1341
1342	END SUBROUTINE routing_simple_initialize
1343
1344
1345	!! ================================================================================================================================
1346	!! SUBROUTINE : routing_simple_main
1347	!!
1348	!>\BRIEF This module routes the water over the continents (runoff and
1349	!! drainage produced by the hydrolc or hydrol module) into the oceans.
1350	!!
1351	!! DESCRIPTION (definitions, functional, design, flags):
1352	!! The routing scheme (Polcher, 2003) carries the water from the runoff and drainage simulated by SECHIBA
1353	!! to the ocean through reservoirs, with some delay. The routing scheme is based on
1354	!! a parametrization of the water flow on a global scale (Miller et al., 1994; Hagemann
1355	!! and Dumenil, 1998). Given the global map of the main watersheds (Oki et al., 1999;
1356	!! Fekete et al., 1999; Vorosmarty et al., 2000) which delineates the boundaries of subbasins
1357	!! and gives the eight possible directions of water flow within the pixel, the surface
1358	!! runoff and the deep drainage are routed to the ocean. The time-step of the routing is one day.
1359	!! The scheme also diagnoses how much water is retained in the foodplains and thus return to soil
1360	!! moisture or is taken out of the rivers for irrigation. \n
1361	!!
1362	!! RECENT CHANGE(S): None
1363	!!
1364	!! MAIN OUTPUT VARIABLE(S):
1365	!! The result of the routing are 3 fluxes :
1366	!! - riverflow : The water which flows out from the major rivers. The flux will be located
1367	!! on the continental grid but this should be a coastal point.
1368	!! - coastalflow : This is the water which flows in a disperse way into the ocean. Essentially these
1369	!! are the outflows from all of the small rivers.
1370	!! - returnflow : This is the water which flows into a land-point - typically rivers which end in
1371	!! the desert. This water will go back into the hydrol module to allow re-evaporation.
1372	!! - irrigation : This is water taken from the reservoir and is being put into the upper
1373	!! layers of the soil.
1374	!! The two first fluxes are in kg/dt and the last two fluxes are in kg/(m^2dt).\n
1375	!!
1376	!! REFERENCE(S) :
1377	!! - Miller JR, Russell GL, Caliri G (1994)
1378	!! Continental-scale river flow in climate models.
1379	!! J. Clim., 7:914-928
1380	!! - Hagemann S and Dumenil L. (1998)
1381	!! A parametrization of the lateral waterflow for the global scale.
1382	!! Clim. Dyn., 14:17-31
1383	!! - Oki, T., T. Nishimura, and P. Dirmeyer (1999)
1384	!! Assessment of annual runoff from land surface models using total runoff integrating pathways (TRIP)
1385	!! J. Meteorol. Soc. Jpn., 77, 235-255
1386	!! - Fekete BM, Charles V, Grabs W (2000)
1387	!! Global, composite runoff fields based on observed river discharge and simulated water balances.
1388	!! Technical report, UNH/GRDC, Global Runoff Data Centre, Koblenz
1389	!! - Vorosmarty, C., B. Fekete, B. Meybeck, and R. Lammers (2000)
1390	!! Global system of rivers: Its role in organizing continental land mass and defining land-to-ocean linkages
1391	!! Global Biogeochem. Cycles, 14, 599-621
1392	!! - Vivant, A-C. (?? 2002)
1393	!! DÃ©veloppement du schÃ©ma de routage et des plaines d'inondation, MSc Thesis, Paris VI University
1394	!! - J. Polcher (2003)
1395	!! Les processus de surface a l'echelle globale et leurs interactions avec l'atmosphere
1396	!! Habilitation a diriger les recherches, Paris VI University, 67pp.
1397	!!
1398	!! FLOWCHART :
1399	!! \latexonly
1400	!! \includegraphics[scale=0.75]{routing_main_flowchart.png}
1401	!! \endlatexonly
1402	!! \n
1403	!_ ================================================================================================================================
1404
1405
1406	SUBROUTINE routing_simple_main(kjit, nbpt, index, &
1407	lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
1408	drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
1409	stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, &
1410	rest_id, hist_id, hist2_id)
1411
1412	IMPLICIT NONE
1413
1414	!! 0 Variable and parameter description
1415	!! 0.1 Input variables
1416
1417	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1418	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
1419	INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless)
1420	INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless)
1421	INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless)
1422	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless)
1423	REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !)
1424	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
1425	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m)
1426	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
1427	REAL(r_std), INTENT(in) :: totfrac_nobio(nbpt) !! Total fraction of no-vegetation (continental ice, lakes ...) (unitless;0-1)
1428	REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) !! Maximal fraction of vegetation (unitless;0-1)
1429	REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt)
1430	REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff (kg/m^2/dt)
1431	REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt)
1432	REAL(r_std), INTENT(in) :: transpot(nbpt,nvm) !! Potential transpiration of the vegetation (kg/m^2/dt)
1433	REAL(r_std), INTENT(in) :: precip_rain(nbpt) !! Rainfall (kg/m^2/dt)
1434	REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
1435	REAL(r_std), INTENT(in) :: humrel(nbpt,nvm) !! Soil moisture stress, root extraction potential (unitless)
1436	REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
1437	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)
1438
1439	!! 0.2 Output variables
1440	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box.
1441	!! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
1442	REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
1443	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)
1444	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)
1445	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)
1446	REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
1447	REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
1448
1449	!! 0.3 Local variables
1450	CHARACTER(LEN=30) :: var_name !! To store variables names for I/O (unitless)
1451	REAL(r_std), DIMENSION(1) :: tmp_day !!
1452	REAL(r_std), DIMENSION(nbpt) :: return_lakes !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
1453	INTEGER(i_std) :: ig, jv !! Indices (unitless)
1454	REAL(r_std), DIMENSION(nbpt) :: tot_vegfrac_nowoody !! Total fraction occupied by grass (0-1,unitless)
1455	REAL(r_std),PARAMETER :: delta_lon=(360./144.)/2.
1456	REAL(r_std),PARAMETER :: delta_lat=(180./144.)/2.
1457
1458	!_ ================================================================================================================================
1459
1460	!
1461	!! Computes the variables averaged between routing time steps and which will be used in subsequent calculations
1462	!
1463	floodout_mean(:) = floodout_mean(:) + floodout(:)
1464	runoff_mean(:) = runoff_mean(:) + runoff(:)
1465	drainage_mean(:) = drainage_mean(:) + drainage(:)
1466	floodtemp(:) = stempdiag(:,floodtemp_lev)
1467	precip_mean(:) = precip_mean(:) + precip_rain(:)
1468	!
1469	!! Computes the total fraction occupied by the grasses and the crops for each grid cell
1470	tot_vegfrac_nowoody(:) = zero
1471	DO jv = 1, nvm
1472	IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
1473	tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv)
1474	END IF
1475	END DO
1476
1477	DO ig = 1, nbpt
1478	IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN
1479	DO jv = 1,nvm
1480	IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
1481	transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig)
1482	END IF
1483	END DO
1484	ELSE
1485	IF (MAXVAL(veget_max(ig,2:nvm)) .GT. min_sechiba) THEN
1486	DO jv = 2, nvm
1487	transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,2:nvm))
1488	ENDDO
1489	ENDIF
1490	ENDIF
1491	ENDDO
1492
1493	!
1494	! Averaged variables (i.e. *dt_sechiba/dt_routing). This accounts for the difference between the shorter
1495	! timestep dt_sechiba of other parts of the model and the long dt_routing timestep (set to one day at present)
1496	!
1497	totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dt_sechiba/dt_routing
1498	k_litt_mean(:) = k_litt_mean(:) + k_litt(:)*dt_sechiba/dt_routing
1499	!
1500	! Only potentially vegetated surfaces are taken into account. At the start of
1501	! the growing seasons we will give more weight to these areas.
1502	!
1503	DO jv=2,nvm
1504	DO ig=1,nbpt
1505	humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)veget_max(ig,jv)dt_sechiba/dt_routing
1506	vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
1507	ENDDO
1508	ENDDO
1509	!
1510	time_counter = time_counter + dt_sechiba
1511	!
1512	! If the time has come we do the routing.
1513	!
1514	IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN
1515	!
1516	! Check the water balance if needed
1517	!
1518	! IF ( check_waterbal ) THEN
1519	! CALL routing_waterbal(nbpt, .TRUE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, &
1520	! & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean)
1521	! ENDIF
1522	!
1523	! Make sure we do not flood north of 49N as there freezing processes start to play a role and they
1524	! are not yet well treated in ORCHIDEE.
1525	!
1526	DO ig=1,nbpt
1527	IF ( lalo(ig,1) > 49.0 ) THEN
1528	floodtemp(ig) = tp_00 - un
1529	ENDIF
1530	ENDDO
1531	!
1532	!! Computes the transport of water in the various reservoirs
1533	!
1534	!ym runoff_mean=0
1535	!ym drainage_mean=0
1536	!ym DO ig=1,nbpt
1537	!ym IF ( lalo(ig,1)-delta_lat < 0 .AND. lalo(ig,1)+delta_lat > 0 .AND. lalo(ig,2)-delta_lon < 32. .AND. lalo(ig,2)+delta_lon > 32.) runoff_mean(ig)=one_day/dt_routing
1538	!ym ENDDO
1539
1540	CALL routing_simple_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, &
1541	& vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, reinf_slope, &
1542	& lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrigation_mean, riverflow_mean, &
1543	& coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res)
1544	!
1545	!! Responsible for storing the water in lakes
1546	!
1547	CALL routing_simple_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, contfrac, return_lakes)
1548	!
1549	returnflow_mean(:) = returnflow_mean(:) + return_lakes(:)
1550	!
1551	!! Check the water balance in the routing scheme
1552	!
1553	! IF ( check_waterbal ) THEN
1554	! CALL routing_waterbal(nbpt, .FALSE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, &
1555	! & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean)
1556	! ENDIF
1557	!
1558	time_counter = zero
1559	!
1560	floodout_mean(:) = zero
1561	runoff_mean(:) = zero
1562	drainage_mean(:) = zero
1563	transpot_mean(:) = zero
1564	precip_mean(:) = zero
1565	!
1566	humrel_mean(:) = zero
1567	totnobio_mean(:) = zero
1568	k_litt_mean(:) = zero
1569	vegtot_mean(:) = zero
1570	!
1571	! Change the units of the routing fluxes from kg/dt_routing into kg/dt_sechiba
1572	! and from m^3/dt_routing into m^3/dt_sechiba
1573	!
1574	returnflow_mean(:) = returnflow_mean(:)/dt_routing*dt_sechiba
1575	reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba
1576	irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba
1577	irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba
1578	!
1579	! Change units as above but at the same time transform the kg/dt_sechiba to m^3/dt_sechiba
1580	!
1581	riverflow_mean(:) = riverflow_mean(:)/dt_routing*dt_sechiba/mille
1582	coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dt_sechiba/mille
1583	hydrographs(:) = hydrographs(:)/dt_routing*dt_sechiba/mille
1584	slowflow_diag(:) = slowflow_diag(:)/dt_routing*dt_sechiba/mille
1585	ENDIF
1586
1587	!
1588	! Return the fraction of routed water for this time step.
1589	!
1590	returnflow(:) = returnflow_mean(:)
1591	reinfiltration(:) = reinfiltration_mean(:)
1592	irrigation(:) = irrigation_mean(:)
1593	riverflow(:) = riverflow_mean(:)
1594	coastalflow(:) = coastalflow_mean(:)
1595
1596	!
1597	! Write diagnostics
1598	!
1599
1600	CALL xios_orchidee_send_field("riversret",returnflow*one_day/dt_sechiba)
1601	CALL xios_orchidee_send_field("hydrographs",hydrographs/dt_sechiba)
1602	CALL xios_orchidee_send_field("slowflow",slowflow_diag/dt_sechiba) ! Qb in m3/s
1603	CALL xios_orchidee_send_field("reinfiltration",reinfiltration)
1604	CALL xios_orchidee_send_field("fastr",fast_diag)
1605	CALL xios_orchidee_send_field("slowr",slow_diag)
1606	CALL xios_orchidee_send_field("streamr",stream_diag)
1607	CALL xios_orchidee_send_field("laker",lake_diag)
1608	CALL xios_orchidee_send_field("pondr",pond_diag)
1609	CALL xios_orchidee_send_field("irrigation",irrigation*one_day/dt_sechiba)
1610	CALL xios_orchidee_send_field("netirrig",irrig_netereq*one_day/dt_sechiba)
1611	CALL xios_orchidee_send_field("floodh",flood_height)
1612	CALL xios_orchidee_send_field("floodr",flood_diag)
1613	delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag
1614	CALL xios_orchidee_send_field("SurfStor",flood_diag+pond_diag+lake_diag)
1615
1616	! Transform from kg/dt_sechiba into m^3/s
1617	CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba)
1618	CALL xios_orchidee_send_field("slowflow",slowflow_diag/mille/dt_sechiba) ! previous id name: Qb
1619	CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
1620	CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
1621
1622
1623	IF ( .NOT. almaoutput ) THEN
1624	!
1625	CALL histwrite_p(hist_id, 'riversret', kjit, returnflow, nbpt, index)
1626	IF (do_floodplains .OR. doponds) THEN
1627	CALL histwrite_p(hist_id, 'reinfiltration', kjit, reinfiltration, nbpt, index)
1628	ENDIF
1629	CALL histwrite_p(hist_id, 'hydrographs', kjit, hydrographs, nbpt, index)
1630	!
1631	CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index)
1632	CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index)
1633	CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index)
1634	IF ( do_floodplains ) THEN
1635	CALL histwrite_p(hist_id, 'floodr', kjit, flood_diag, nbpt, index)
1636	CALL histwrite_p(hist_id, 'floodh', kjit, flood_height, nbpt, index)
1637	ENDIF
1638	CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, nbpt, index)
1639	CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index)
1640	!
1641	IF ( do_irrigation ) THEN
1642	CALL histwrite_p(hist_id, 'irrigation', kjit, irrigation, nbpt, index)
1643	CALL histwrite_p(hist_id, 'returnflow', kjit, returnflow, nbpt, index)
1644	CALL histwrite_p(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
1645	ENDIF
1646	!
1647	ELSE
1648	!
1649	delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag
1650	CALL histwrite_p(hist_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index)
1651	CALL histwrite_p(hist_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index)
1652	CALL histwrite_p(hist_id, 'Dis', kjit, hydrographs, nbpt, index)
1653	IF ( do_irrigation ) THEN
1654	CALL histwrite_p(hist_id, 'Qirrig', kjit, irrigation, nbpt, index)
1655	CALL histwrite_p(hist_id, 'Qirrig_req', kjit, irrig_netereq, nbpt, index)
1656	ENDIF
1657	!
1658	ENDIF
1659	IF ( hist2_id > 0 ) THEN
1660	IF ( .NOT. almaoutput ) THEN
1661	!
1662	CALL histwrite_p(hist2_id, 'riversret', kjit, returnflow, nbpt, index)
1663	IF (do_floodplains .OR. doponds) THEN
1664	CALL histwrite_p(hist2_id, 'reinfiltration', kjit, reinfiltration, nbpt, index)
1665	ENDIF
1666	CALL histwrite_p(hist2_id, 'hydrographs', kjit, hydrographs, nbpt, index)
1667	!
1668	CALL histwrite_p(hist2_id, 'fastr', kjit, fast_diag, nbpt, index)
1669	CALL histwrite_p(hist2_id, 'slowr', kjit, slow_diag, nbpt, index)
1670	IF ( do_floodplains ) THEN
1671	CALL histwrite_p(hist2_id, 'floodr', kjit, flood_diag, nbpt, index)
1672	CALL histwrite_p(hist2_id, 'floodh', kjit, flood_height, nbpt, index)
1673	ENDIF
1674	CALL histwrite_p(hist2_id, 'pondr', kjit, pond_diag, nbpt, index)
1675	CALL histwrite_p(hist2_id, 'streamr', kjit, stream_diag, nbpt, index)
1676	CALL histwrite_p(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index)
1677	!
1678	IF ( do_irrigation ) THEN
1679	CALL histwrite_p(hist2_id, 'irrigation', kjit, irrigation, nbpt, index)
1680	CALL histwrite_p(hist2_id, 'returnflow', kjit, returnflow, nbpt, index)
1681	CALL histwrite_p(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
1682	ENDIF
1683	!
1684	ELSE
1685	!
1686	delsurfstor=delsurfstor + flood_diag + pond_diag + lake_diag
1687	CALL histwrite_p(hist2_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index)
1688	CALL histwrite_p(hist2_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index)
1689	CALL histwrite_p(hist2_id, 'Dis', kjit, hydrographs, nbpt, index)
1690	!
1691	ENDIF
1692	ENDIF
1693	!
1694	!
1695	END SUBROUTINE routing_simple_main
1696
1697	!! =============================================================================================================================
1698	!! SUBROUTINE: routing_simple_finalize
1699	!!
1700	!>\BRIEF Write to restart file
1701	!!
1702	!! DESCRIPTION: Write module variables to restart file
1703	!!
1704	!! RECENT CHANGE(S)
1705	!!
1706	!! REFERENCE(S)
1707	!!
1708	!! FLOWCHART
1709	!! \n
1710	!_ ==============================================================================================================================
1711
1712	SUBROUTINE routing_simple_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
1713	USE xios
1714	IMPLICIT NONE
1715
1716	!! 0 Variable and parameter description
1717	!! 0.1 Input variables
1718	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1719	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
1720	INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless)
1721	REAL(r_std), INTENT(in) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1)
1722	REAL(r_std), INTENT(in) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg)
1723
1724	!! 0.2 Local variables
1725	REAL(r_std), DIMENSION(1) :: tmp_day
1726
1727	!_ ================================================================================================================================
1728
1729	IF (is_omp_root) THEN
1730	CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r)
1731	CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r)
1732	CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r)
1733	ENDIF
1734
1735	END SUBROUTINE routing_simple_finalize
1736
1737	!! ================================================================================================================================
1738	!! SUBROUTINE : routing_simple_init_1
1739	!!
1740	!>\BRIEF This subroutine allocates the memory and get the fixed fields from the restart file.
1741	!!
1742	!! DESCRIPTION: Privat subroutine to the module routing_simple. This subroutine is called in the begining
1743	!! of routing_simple_initialize
1744	!!
1745	!! RECENT CHANGE(S): None
1746	!!
1747	!! MAIN OUTPUT VARIABLE(S):
1748	!!
1749	!! REFERENCES : None
1750	!!
1751	!! FLOWCHART :None
1752	!! \n
1753	!_ ================================================================================================================================
1754
1755	SUBROUTINE routing_simple_init_1(kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
1756	riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id)
1757
1758	IMPLICIT NONE
1759
1760	!! 0 Variable and parameter description
1761	!! 0.1 Input variables
1762
1763	INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless)
1764	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
1765	INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index !! Indices of the points on the map (unitless)
1766	REAL(r_std), DIMENSION(nbpt,nslm),INTENT(in) :: stempdiag !! Temperature profile in soil
1767	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless)
1768
1769	!! 0.2 Output variables
1770	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: returnflow !! The water flow from lakes and swamps which returns into the grid box.
1771	!! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
1772	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: reinfiltration !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
1773	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)
1774	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)
1775	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)
1776	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_frac !! Flooded fraction of the grid box (unitless;0-1)
1777	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_res !! Diagnostic of water amount in the floodplains reservoir (kg)
1778
1779	!! 0.3 Local variables
1780	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
1781	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g !! A temporary real array for the integers
1782	REAL(r_std), DIMENSION(1) :: tmp_day !!
1783	REAL(r_std) :: ratio !! Diagnostic ratio to check that dt_routing is a multiple of dt_sechiba (unitless)
1784	REAL(r_std) :: totarea !! Total area of basin (m^2)
1785	INTEGER(i_std) :: ier, ig, ib, ipn(1) !! Indices (unitless)
1786
1787	!_ ================================================================================================================================
1788	!
1789	!
1790	! These variables will require the configuration infrastructure
1791	!
1792	!Config Key = DT_ROUTING
1793	!Config If = RIVER_ROUTING
1794	!Config Desc = Time step of the routing scheme
1795	!Config Def = one_day
1796	!Config Help = This values gives the time step in seconds of the routing scheme.
1797	!Config It should be multiple of the main time step of ORCHIDEE. One day
1798	!Config is a good value.
1799	!Config Units = [seconds]
1800	!
1801	dt_routing = one_day
1802	CALL getin_p('DT_ROUTING', dt_routing)
1803	!
1804	!Config Key = ROUTING_RIVERS
1805	!Config If = RIVER_ROUTING
1806	!Config Desc = Number of rivers
1807	!Config Def = 50
1808	!Config Help = This parameter chooses the number of largest river basins
1809	!Config which should be treated as independently as rivers and not
1810	!Config flow into the oceans as diffusion coastal flow.
1811	!Config Units = [-]
1812	num_largest = 50
1813	CALL getin_p('ROUTING_RIVERS', num_largest)
1814	!
1815	!Config Key = DO_FLOODINFILT
1816	!Config Desc = Should floodplains reinfiltrate into the soil
1817	!Config If = RIVER_ROUTING
1818	!Config Def = n
1819	!Config Help = This parameters allows the user to ask the model
1820	!Config to take into account the flood plains reinfiltration
1821	!Config into the soil moisture. It then can go
1822	!Config back to the slow and fast reservoirs
1823	!Config Units = [FLAG]
1824	!
1825	dofloodinfilt = .FALSE.
1826	CALL getin_p('DO_FLOODINFILT', dofloodinfilt)
1827	!
1828	!Config Key = DO_SWAMPS
1829	!Config Desc = Should we include swamp parameterization
1830	!Config If = RIVER_ROUTING
1831	!Config Def = n
1832	!Config Help = This parameters allows the user to ask the model
1833	!Config to take into account the swamps and return
1834	!Config the water into the bottom of the soil. It then can go
1835	!Config back to the atmopshere. This tried to simulate
1836	!Config internal deltas of rivers.
1837	!Config Units = [FLAG]
1838	!
1839	doswamps = .FALSE.
1840	CALL getin_p('DO_SWAMPS', doswamps)
1841	!
1842	!Config Key = DO_PONDS
1843	!Config Desc = Should we include ponds
1844	!Config If = RIVER_ROUTING
1845	!Config Def = n
1846	!Config Help = This parameters allows the user to ask the model
1847	!Config to take into account the ponds and return
1848	!Config the water into the soil moisture. It then can go
1849	!Config back to the atmopshere. This tried to simulate
1850	!Config little ponds especially in West Africa.
1851	!Config Units = [FLAG]
1852	!
1853	doponds = .FALSE.
1854	CALL getin_p('DO_PONDS', doponds)
1855
1856	!Config Key = SLOW_TCST
1857	!Config Desc = Time constant for the slow reservoir
1858	!Config If = RIVER_ROUTING
1859	!Config Def = 25.0
1860	!Config Help = This parameters allows the user to fix the
1861	!Config time constant (in days) of the slow reservoir
1862	!Config in order to get better river flows for
1863	!Config particular regions.
1864	!Config Units = [days]
1865	!
1866	!> A value for property of each reservoir (in day/m) is given to compute a time constant (in day)
1867	!> for each reservoir (product of tcst and topo_resid).
1868	!> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only,
1869	!> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and
1870	!> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir"
1871	!> have the highest value in order to simulate the groundwater.
1872	!> The "stream reservoir", which represents all the water of the stream, has the lowest value.
1873	!> Those figures are the same for all the basins of the world.
1874	CALL getin_p('SLOW_TCST', slow_tcst)
1875	!
1876	!Config Key = FAST_TCST
1877	!Config Desc = Time constant for the fast reservoir
1878	!Config If = RIVER_ROUTING
1879	!Config Def = 3.0
1880	!Config Help = This parameters allows the user to fix the
1881	!Config time constant (in days) of the fast reservoir
1882	!Config in order to get better river flows for
1883	!Config particular regions.
1884	!Config Units = [days]
1885	CALL getin_p('FAST_TCST', fast_tcst)
1886	!
1887	!Config Key = STREAM_TCST
1888	!Config Desc = Time constant for the stream reservoir
1889	!Config If = RIVER_ROUTING
1890	!Config Def = 0.24
1891	!Config Help = This parameters allows the user to fix the
1892	!Config time constant (in days) of the stream reservoir
1893	!Config in order to get better river flows for
1894	!Config particular regions.
1895	!Config Units = [days]
1896	CALL getin_p('STREAM_TCST', stream_tcst)
1897	!
1898	!Config Key = FLOOD_TCST
1899	!Config Desc = Time constant for the flood reservoir
1900	!Config If = RIVER_ROUTING
1901	!Config Def = 4.0
1902	!Config Help = This parameters allows the user to fix the
1903	!Config time constant (in days) of the flood reservoir
1904	!Config in order to get better river flows for
1905	!Config particular regions.
1906	!Config Units = [days]
1907	CALL getin_p('FLOOD_TCST', flood_tcst)
1908
1909	!Config Key = SWAMP_CST
1910	!Config Desc = Fraction of the river that flows back to swamps
1911	!Config If = RIVER_ROUTING
1912	!Config Def = 0.2
1913	!Config Help = This parameters allows the user to fix the
1914	!Config fraction of the river transport
1915	!Config that flows to swamps
1916	!Config Units = [-]
1917	CALL getin_p('SWAMP_CST', swamp_cst)
1918
1919	!Config Key = FLOOD_BETA
1920	!Config Desc = Parameter to fix the shape of the floodplain
1921	!Config If = RIVER_ROUTING
1922	!Config Def = 2.0
1923	!Config Help = Parameter to fix the shape of the floodplain
1924	!Config (>1 for convex edges, <1 for concave edges)
1925	!Config Units = [-]
1926	CALL getin_p("FLOOD_BETA", beta)
1927
1928	!Config Key = POND_BETAP
1929	!Config Desc = Ratio of the basin surface intercepted by ponds and the maximum surface of ponds
1930	!Config If = RIVER_ROUTING
1931	!Config Def = 0.5
1932	!Config Help =
1933	!Config Units = [-]
1934	CALL getin_p("POND_BETAP", betap)
1935
1936	!Config Key = FLOOD_CRI
1937	!Config Desc = Potential height for which all the basin is flooded
1938	!Config If = DO_FLOODPLAINS or DO_PONDS
1939	!Config Def = 2000.
1940	!Config Help =
1941	!Config Units = [mm]
1942	CALL getin_p("FLOOD_CRI", floodcri)
1943
1944	!Config Key = POND_CRI
1945	!Config Desc = Potential height for which all the basin is a pond
1946	!Config If = DO_FLOODPLAINS or DO_PONDS
1947	!Config Def = 2000.
1948	!Config Help =
1949	!Config Units = [mm]
1950	CALL getin_p("POND_CRI", pondcri)
1951
1952
1953	! In order to simplify the time cascade check that dt_routing
1954	! is a multiple of dt_sechiba
1955	!
1956	ratio = dt_routing/dt_sechiba
1957	IF ( ABS(NINT(ratio) - ratio) .GT. 10*EPSILON(ratio)) THEN
1958	WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
1959	WRITE(numout,*) "The chosen time step for the routing is not a multiple of the"
1960	WRITE(numout,*) "main time step of the model. We will change dt_routing so that"
1961	WRITE(numout,*) "this condition os fulfilled"
1962	dt_routing = NINT(ratio) * dt_sechiba
1963	WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
1964	ENDIF
1965	!
1966	IF ( dt_routing .LT. dt_sechiba) THEN
1967	WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
1968	WRITE(numout,*) 'The routing timestep can not be smaller than the one'
1969	WRITE(numout,*) 'of the model. We reset its value to the model''s timestep.'
1970	WRITE(numout,*) 'The old DT_ROUTING is : ', dt_routing
1971	dt_routing = dt_sechiba
1972	WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
1973	ENDIF
1974	!
1975	var_name ="routingcounter"
1976	IF (is_root_prc) THEN
1977	CALL ioconf_setatt('UNITS', 's')
1978	CALL ioconf_setatt('LONG_NAME','Time counter for the routing scheme')
1979	CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day)
1980	IF (tmp_day(1) == val_exp) THEN
1981	time_counter = zero
1982	ELSE
1983	time_counter = tmp_day(1)
1984	ENDIF
1985	ENDIF
1986	CALL bcast(time_counter)
1987
1988
1989	ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier)
1990	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_loc','','')
1991
1992	ALLOCATE (routing_area_glo(nbp_glo,nbasmax))
1993	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_glo','','')
1994	var_name = 'routingarea'
1995	IF (is_root_prc) THEN
1996	CALL ioconf_setatt('UNITS', 'm^2')
1997	CALL ioconf_setatt('LONG_NAME','Area of basin')
1998	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g)
1999	ENDIF
2000	CALL scatter(routing_area_glo,routing_area_loc)
2001	routing_area=>routing_area_loc
2002
2003	ALLOCATE (tmp_real_g(nbp_glo,nbasmax), stat=ier)
2004	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for tmp_real_g','','')
2005
2006	ALLOCATE (route_togrid_loc(nbpt,nbasmax), stat=ier)
2007	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_loc','','')
2008	ALLOCATE (route_togrid_glo(nbp_glo,nbasmax), stat=ier) ! used in global in routing_simple_flow
2009	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_glo','','')
2010
2011	IF (is_root_prc) THEN
2012	var_name = 'routetogrid'
2013	CALL ioconf_setatt('UNITS', '-')
2014	CALL ioconf_setatt('LONG_NAME','Grid into which the basin flows')
2015	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
2016	route_togrid_glo(:,:) = undef_int
2017	WHERE ( tmp_real_g .LT. val_exp )
2018	route_togrid_glo = NINT(tmp_real_g)
2019	ENDWHERE
2020	ENDIF
2021	CALL bcast(route_togrid_glo) ! used in global in routing_simple_flow
2022	CALL scatter(route_togrid_glo,route_togrid_loc)
2023	route_togrid=>route_togrid_loc
2024	!
2025	ALLOCATE (route_tobasin_loc(nbpt,nbasmax), stat=ier)
2026	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_loc','','')
2027
2028	ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax), stat=ier)
2029	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_glo','','')
2030
2031	IF (is_root_prc) THEN
2032	var_name = 'routetobasin'
2033	CALL ioconf_setatt('UNITS', '-')
2034	CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes')
2035	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
2036	route_tobasin_glo = undef_int
2037	WHERE ( tmp_real_g .LT. val_exp )
2038	route_tobasin_glo = NINT(tmp_real_g)
2039	ENDWHERE
2040	ENDIF
2041	CALL scatter(route_tobasin_glo,route_tobasin_loc)
2042	route_tobasin=>route_tobasin_loc
2043	!
2044	! nbintobasin
2045	!
2046	ALLOCATE (route_nbintobas_loc(nbpt,nbasmax), stat=ier)
2047	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_loc','','')
2048	ALLOCATE (route_nbintobas_glo(nbp_glo,nbasmax), stat=ier)
2049	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_glo','','')
2050
2051	IF (is_root_prc) THEN
2052	var_name = 'routenbintobas'
2053	CALL ioconf_setatt('UNITS', '-')
2054	CALL ioconf_setatt('LONG_NAME','Number of basin into current one')
2055	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
2056	route_nbintobas_glo = undef_int
2057	WHERE ( tmp_real_g .LT. val_exp )
2058	route_nbintobas_glo = NINT(tmp_real_g)
2059	ENDWHERE
2060	ENDIF
2061	CALL scatter(route_nbintobas_glo,route_nbintobas_loc)
2062	route_nbintobas=>route_nbintobas_loc
2063	!
2064	ALLOCATE (global_basinid_loc(nbpt,nbasmax), stat=ier)
2065	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_loc','','')
2066	ALLOCATE (global_basinid_glo(nbp_glo,nbasmax), stat=ier)
2067	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_glo','','')
2068
2069	IF (is_root_prc) THEN
2070	var_name = 'basinid'
2071	CALL ioconf_setatt('UNITS', '-')
2072	CALL ioconf_setatt('LONG_NAME','ID of basin')
2073	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
2074	global_basinid_glo = undef_int
2075	WHERE ( tmp_real_g .LT. val_exp )
2076	global_basinid_glo = NINT(tmp_real_g)
2077	ENDWHERE
2078	ENDIF
2079	CALL scatter(global_basinid_glo,global_basinid_loc)
2080	global_basinid=>global_basinid_loc
2081	!
2082	ALLOCATE (topo_resid_loc(nbpt,nbasmax), stat=ier)
2083	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_loc','','')
2084	ALLOCATE (topo_resid_glo(nbp_glo,nbasmax), stat=ier)
2085	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_glo','','')
2086
2087	IF (is_root_prc) THEN
2088	var_name = 'topoindex'
2089	CALL ioconf_setatt('UNITS', 'm')
2090	CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time')
2091	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g)
2092	ENDIF
2093	CALL scatter(topo_resid_glo,topo_resid_loc)
2094	topo_resid=>topo_resid_loc
2095
2096	ALLOCATE (fast_reservoir(nbpt,nbasmax), stat=ier)
2097	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_reservoir','','')
2098	var_name = 'fastres'
2099	CALL ioconf_setatt_p('UNITS', 'Kg')
2100	CALL ioconf_setatt_p('LONG_NAME','Water in the fast reservoir')
2101	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g)
2102	CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero)
2103
2104	ALLOCATE (slow_reservoir(nbpt,nbasmax), stat=ier)
2105	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slow_reservoir','','')
2106	var_name = 'slowres'
2107	CALL ioconf_setatt_p('UNITS', 'Kg')
2108	CALL ioconf_setatt_p('LONG_NAME','Water in the slow reservoir')
2109	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g)
2110	CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero)
2111
2112	ALLOCATE (stream_reservoir(nbpt,nbasmax), stat=ier)
2113	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for stream_reservoir','','')
2114	var_name = 'streamres'
2115	CALL ioconf_setatt_p('UNITS', 'Kg')
2116	CALL ioconf_setatt_p('LONG_NAME','Water in the stream reservoir')
2117	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g)
2118	CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero)
2119
2120	ALLOCATE (flood_reservoir(nbpt,nbasmax), stat=ier)
2121	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_reservoir','','')
2122	var_name = 'floodres'
2123	CALL ioconf_setatt_p('UNITS', 'Kg')
2124	CALL ioconf_setatt_p('LONG_NAME','Water in the flood reservoir')
2125	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_reservoir, "gather", nbp_glo, index_g)
2126	CALL setvar_p (flood_reservoir, val_exp, 'NO_KEYWORD', zero)
2127
2128	ALLOCATE (flood_frac_bas(nbpt,nbasmax), stat=ier)
2129	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_frac_bas','','')
2130	var_name = 'flood_frac_bas'
2131	CALL ioconf_setatt_p('UNITS', '-')
2132	CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per basin')
2133	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_frac_bas, "gather", nbp_glo, index_g)
2134	CALL setvar_p (flood_frac_bas, val_exp, 'NO_KEYWORD', zero)
2135
2136	ALLOCATE (flood_height(nbpt), stat=ier)
2137	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_height','','')
2138	var_name = 'floodh'
2139	CALL ioconf_setatt_p('UNITS', '-')
2140	CALL ioconf_setatt_p('LONG_NAME','')
2141	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g)
2142	CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero)
2143
2144	ALLOCATE (pond_frac(nbpt), stat=ier)
2145	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for pond_frac','','')
2146	var_name = 'pond_frac'
2147	CALL ioconf_setatt_p('UNITS', '-')
2148	CALL ioconf_setatt_p('LONG_NAME','Pond fraction per grid box')
2149	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g)
2150	CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero)
2151
2152	var_name = 'flood_frac'
2153	CALL ioconf_setatt_p('UNITS', '-')
2154	CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per grid box')
2155	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g)
2156	CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero)
2157
2158	var_name = 'flood_res'
2159	CALL ioconf_setatt_p('UNITS','mm')
2160	CALL ioconf_setatt_p('LONG_NAME','Flooded quantity (estimation)')
2161	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_res, "gather", nbp_glo, index_g)
2162	CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero)
2163	! flood_res = zero
2164
2165	ALLOCATE (lake_reservoir(nbpt), stat=ier)
2166	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lake_reservoir','','')
2167	var_name = 'lakeres'
2168	CALL ioconf_setatt_p('UNITS', 'Kg')
2169	CALL ioconf_setatt_p('LONG_NAME','Water in the lake reservoir')
2170	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g)
2171	CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero)
2172
2173	ALLOCATE (pond_reservoir(nbpt), stat=ier)
2174	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for pond_reservoir','','')
2175	var_name = 'pondres'
2176	CALL ioconf_setatt_p('UNITS', 'Kg')
2177	CALL ioconf_setatt_p('LONG_NAME','Water in the pond reservoir')
2178	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_reservoir, "gather", nbp_glo, index_g)
2179	CALL setvar_p (pond_reservoir, val_exp, 'NO_KEYWORD', zero)
2180	!
2181	! Map of irrigated areas
2182	!
2183	IF ( do_irrigation ) THEN
2184	ALLOCATE (irrigated(nbpt), stat=ier)
2185	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrigated','','')
2186	var_name = 'irrigated'
2187	CALL ioconf_setatt_p('UNITS', 'm^2')
2188	CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area')
2189	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g)
2190	CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba)
2191	ENDIF
2192
2193	IF ( do_floodplains ) THEN
2194	ALLOCATE (floodplains(nbpt), stat=ier)
2195	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodplains','','')
2196	var_name = 'floodplains'
2197	CALL ioconf_setatt_p('UNITS', 'm^2')
2198	CALL ioconf_setatt_p('LONG_NAME','Surface which can be flooded')
2199	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodplains, "gather", nbp_glo, index_g)
2200	CALL setvar_p (floodplains, val_exp, 'NO_KEYWORD', undef_sechiba)
2201	ENDIF
2202	IF ( doswamps ) THEN
2203	ALLOCATE (swamp(nbpt), stat=ier)
2204	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for swamp','','')
2205	var_name = 'swamp'
2206	CALL ioconf_setatt_p('UNITS', 'm^2')
2207	CALL ioconf_setatt_p('LONG_NAME','Surface which can become swamp')
2208	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., swamp, "gather", nbp_glo, index_g)
2209	CALL setvar_p (swamp, val_exp, 'NO_KEYWORD', undef_sechiba)
2210	ENDIF
2211	!
2212	! Put into the restart file the fluxes so that they can be regenerated at restart.
2213	!
2214	ALLOCATE (lakeinflow_mean(nbpt), stat=ier)
2215	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lakeinflow_mean','','')
2216	var_name = 'lakeinflow'
2217	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
2218	CALL ioconf_setatt_p('LONG_NAME','Lake inflow')
2219	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
2220	CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
2221
2222	ALLOCATE (returnflow_mean(nbpt), stat=ier)
2223	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for returnflow_mean','','')
2224	var_name = 'returnflow'
2225	CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
2226	CALL ioconf_setatt_p('LONG_NAME','Deep return flux')
2227	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g)
2228	CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero)
2229	returnflow(:) = returnflow_mean(:)
2230
2231	ALLOCATE (reinfiltration_mean(nbpt), stat=ier)
2232	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for reinfiltration_mean','','')
2233	var_name = 'reinfiltration'
2234	CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
2235	CALL ioconf_setatt_p('LONG_NAME','Top return flux')
2236	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinfiltration_mean, "gather", nbp_glo, index_g)
2237	CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero)
2238	reinfiltration(:) = reinfiltration_mean(:)
2239
2240	ALLOCATE (irrigation_mean(nbpt), stat=ier)
2241	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrigation_mean','','')
2242	ALLOCATE (irrig_netereq(nbpt), stat=ier)
2243	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrig_netereq','','')
2244	irrig_netereq(:) = zero
2245
2246	IF ( do_irrigation ) THEN
2247	var_name = 'irrigation'
2248	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
2249	CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux')
2250	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g)
2251	CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero)
2252	ELSE
2253	irrigation_mean(:) = zero
2254	ENDIF
2255	irrigation(:) = irrigation_mean(:)
2256
2257	ALLOCATE (riverflow_mean(nbpt), stat=ier)
2258	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for riverflow_mean','','')
2259	var_name = 'riverflow'
2260	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
2261	CALL ioconf_setatt_p('LONG_NAME','River flux into the sea')
2262	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
2263	CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
2264	riverflow(:) = riverflow_mean(:)
2265
2266	ALLOCATE (coastalflow_mean(nbpt), stat=ier)
2267	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for coastalflow_mean','','')
2268	var_name = 'coastalflow'
2269	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
2270	CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea')
2271	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
2272	CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
2273	coastalflow(:) = coastalflow_mean(:)
2274
2275	! Locate it at the 2m level
2276	ipn = MINLOC(ABS(diaglev-2))
2277	floodtemp_lev = ipn(1)
2278	ALLOCATE (floodtemp(nbpt), stat=ier)
2279	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodtemp','','')
2280	floodtemp(:) = stempdiag(:,floodtemp_lev)
2281
2282	ALLOCATE(hydrographs(nbpt), stat=ier)
2283	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrographs','','')
2284	var_name = 'hydrographs'
2285	CALL ioconf_setatt_p('UNITS', 'm^3/dt')
2286	CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid')
2287	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g)
2288	CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero)
2289
2290	ALLOCATE(slowflow_diag(nbpt), stat=ier)
2291	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slowflow_diag','','')
2292	var_name = 'slowflow_diag'
2293	CALL ioconf_setatt_p('UNITS', 'm^3/dt')
2294	CALL ioconf_setatt_p('LONG_NAME','Slowflow hydrograph at outlow of grid')
2295	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,slowflow_diag, "gather", nbp_glo, index_g)
2296	CALL setvar_p (slowflow_diag, val_exp, 'NO_KEYWORD', zero)
2297
2298	!
2299	! The diagnostic variables, they are initialized from the above restart variables.
2300	!
2301	ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), flood_diag(nbpt), &
2302	& pond_diag(nbpt), lake_diag(nbpt), delsurfstor(nbpt), stat=ier)
2303	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_diag,..','','')
2304
2305	fast_diag(:) = zero
2306	slow_diag(:) = zero
2307	stream_diag(:) = zero
2308	flood_diag(:) = zero
2309	pond_diag(:) = zero
2310	lake_diag(:) = zero
2311	delsurfstor(:) = zero
2312
2313	DO ig=1,nbpt
2314	totarea = zero
2315	DO ib=1,nbasmax
2316	totarea = totarea + routing_area(ig,ib)
2317	fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib)
2318	slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib)
2319	stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib)
2320	flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib)
2321	ENDDO
2322	!
2323	fast_diag(ig) = fast_diag(ig)/totarea
2324	slow_diag(ig) = slow_diag(ig)/totarea
2325	stream_diag(ig) = stream_diag(ig)/totarea
2326	flood_diag(ig) = flood_diag(ig)/totarea
2327	!
2328	! This is the volume of the lake scaled to the entire grid.
2329	! It would be batter to scale it to the size of the lake
2330	! but this information is not yet available.
2331	!
2332	lake_diag(ig) = lake_reservoir(ig)/totarea
2333	!
2334	ENDDO
2335	!
2336	! Get from the restart the fluxes we accumulated.
2337	!
2338	ALLOCATE (floodout_mean(nbpt), stat=ier)
2339	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodout_mean','','')
2340	var_name = 'floodout_route'
2341	CALL ioconf_setatt_p('UNITS', 'Kg')
2342	CALL ioconf_setatt_p('LONG_NAME','Accumulated flow out of floodplains for routing')
2343	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g)
2344	CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero)
2345
2346	ALLOCATE (runoff_mean(nbpt), stat=ier)
2347	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for runoff_mean','','')
2348	var_name = 'runoff_route'
2349	CALL ioconf_setatt_p('UNITS', 'Kg')
2350	CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing')
2351	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
2352	CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
2353
2354	ALLOCATE(drainage_mean(nbpt), stat=ier)
2355	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for drainage_mean','','')
2356	var_name = 'drainage_route'
2357	CALL ioconf_setatt_p('UNITS', 'Kg')
2358	CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing')
2359	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
2360	CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
2361
2362	ALLOCATE(transpot_mean(nbpt), stat=ier)
2363	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for transpot_mean','','')
2364	var_name = 'transpot_route'
2365	CALL ioconf_setatt_p('UNITS', 'Kg/m^2')
2366	CALL ioconf_setatt_p('LONG_NAME','Accumulated potential transpiration for routing/irrigation')
2367	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., transpot_mean, "gather", nbp_glo, index_g)
2368	CALL setvar_p (transpot_mean, val_exp, 'NO_KEYWORD', zero)
2369
2370	ALLOCATE(precip_mean(nbpt), stat=ier)
2371	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for precip_mean','','')
2372	var_name = 'precip_route'
2373	CALL ioconf_setatt_p('UNITS', 'Kg/m^2')
2374	CALL ioconf_setatt_p('LONG_NAME','Accumulated rain precipitation for irrigation')
2375	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g)
2376	CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero)
2377
2378	ALLOCATE(humrel_mean(nbpt), stat=ier)
2379	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for humrel_mean','','')
2380	var_name = 'humrel_route'
2381	CALL ioconf_setatt_p('UNITS', '-')
2382	CALL ioconf_setatt_p('LONG_NAME','Mean humrel for irrigation')
2383	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g)
2384	CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un)
2385
2386	ALLOCATE(k_litt_mean(nbpt), stat=ier)
2387	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for k_litt_mean','','')
2388	var_name = 'k_litt_route'
2389	CALL ioconf_setatt_p('UNITS', '-')
2390	CALL ioconf_setatt_p('LONG_NAME','Mean cond. for litter')
2391	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g)
2392	CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero)
2393
2394	ALLOCATE(totnobio_mean(nbpt), stat=ier)
2395	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for totnobio_mean','','')
2396	var_name = 'totnobio_route'
2397	CALL ioconf_setatt_p('UNITS', '-')
2398	CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of no bio for irrigation')
2399	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g)
2400	CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero)
2401
2402	ALLOCATE(vegtot_mean(nbpt), stat=ier)
2403	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for vegtot_mean','','')
2404	var_name = 'vegtot_route'
2405	CALL ioconf_setatt_p('UNITS', '-')
2406	CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of vegetation')
2407	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g)
2408	CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un)
2409	!
2410	!
2411	DEALLOCATE(tmp_real_g)
2412	!
2413	! Allocate diagnostic variables
2414	!
2415	ALLOCATE(hydrodiag_loc(nbpt,nbasmax),hydrodiag_glo(nbp_glo,nbasmax),stat=ier)
2416	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrodiag_glo','','')
2417	hydrodiag=>hydrodiag_loc
2418
2419	ALLOCATE(hydroupbasin_loc(nbpt),hydroupbasin_glo(nbp_glo), stat=ier)
2420	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydroupbasin_glo','','')
2421	hydroupbasin=>hydroupbasin_loc
2422
2423	END SUBROUTINE routing_simple_init_1
2424
2425
2426	!! ================================================================================================================================
2427	!! SUBROUTINE : routing_simple_clear
2428	!!
2429	!>\BRIEF This subroutine deallocates the block memory previously allocated.
2430	!!
2431	!! DESCRIPTION: This subroutine deallocates the block memory previously allocated.
2432	!!
2433	!! RECENT CHANGE(S): None
2434	!!
2435	!! MAIN OUTPUT VARIABLE(S):
2436	!!
2437	!! REFERENCES : None
2438	!!
2439	!! FLOWCHART :None
2440	!! \n
2441	!_ ================================================================================================================================
2442
2443	SUBROUTINE routing_simple_clear()
2444
2445	IF (is_omp_root) THEN
2446	IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r)
2447	IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1)
2448	IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r)
2449	IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r)
2450	IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r)
2451	IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r)
2452	IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r)
2453	ENDIF
2454
2455
2456	IF (ALLOCATED(routing_area_loc)) DEALLOCATE(routing_area_loc)
2457	IF (ALLOCATED(route_togrid_loc)) DEALLOCATE(route_togrid_loc)
2458	IF (ALLOCATED(route_tobasin_loc)) DEALLOCATE(route_tobasin_loc)
2459	IF (ALLOCATED(route_nbintobas_loc)) DEALLOCATE(route_nbintobas_loc)
2460	IF (ALLOCATED(global_basinid_loc)) DEALLOCATE(global_basinid_loc)
2461	IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc)
2462	IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo)
2463	IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo)
2464	IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo)
2465	IF (ALLOCATED(route_nbintobas_glo)) DEALLOCATE(route_nbintobas_glo)
2466	IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo)
2467	IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo)
2468	IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir)
2469	IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir)
2470	IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir)
2471	IF (ALLOCATED(flood_reservoir)) DEALLOCATE(flood_reservoir)
2472	IF (ALLOCATED(flood_frac_bas)) DEALLOCATE(flood_frac_bas)
2473	IF (ALLOCATED(flood_height)) DEALLOCATE(flood_height)
2474	IF (ALLOCATED(pond_frac)) DEALLOCATE(pond_frac)
2475	IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir)
2476	IF (ALLOCATED(pond_reservoir)) DEALLOCATE(pond_reservoir)
2477	IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_mean)
2478	IF (ALLOCATED(reinfiltration_mean)) DEALLOCATE(reinfiltration_mean)
2479	IF (ALLOCATED(riverflow_mean)) DEALLOCATE(riverflow_mean)
2480	IF (ALLOCATED(coastalflow_mean)) DEALLOCATE(coastalflow_mean)
2481	IF (ALLOCATED(lakeinflow_mean)) DEALLOCATE(lakeinflow_mean)
2482	IF (ALLOCATED(runoff_mean)) DEALLOCATE(runoff_mean)
2483	IF (ALLOCATED(floodout_mean)) DEALLOCATE(floodout_mean)
2484	IF (ALLOCATED(drainage_mean)) DEALLOCATE(drainage_mean)
2485	IF (ALLOCATED(transpot_mean)) DEALLOCATE(transpot_mean)
2486	IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean)
2487	IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean)
2488	IF (ALLOCATED(k_litt_mean)) DEALLOCATE(k_litt_mean)
2489	IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean)
2490	IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean)
2491	IF (ALLOCATED(floodtemp)) DEALLOCATE(floodtemp)
2492	IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc)
2493	IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo)
2494	IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc)
2495	IF (ALLOCATED(hydroupbasin_glo)) DEALLOCATE(hydroupbasin_glo)
2496	IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs)
2497	IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag)
2498	IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean)
2499	IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated)
2500	IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains)
2501	IF (ALLOCATED(swamp)) DEALLOCATE(swamp)
2502	IF (ALLOCATED(fast_diag)) DEALLOCATE(fast_diag)
2503	IF (ALLOCATED(slow_diag)) DEALLOCATE(slow_diag)
2504	IF (ALLOCATED(stream_diag)) DEALLOCATE(stream_diag)
2505	IF (ALLOCATED(flood_diag)) DEALLOCATE(flood_diag)
2506	IF (ALLOCATED(pond_diag)) DEALLOCATE(pond_diag)
2507	IF (ALLOCATED(lake_diag)) DEALLOCATE(lake_diag)
2508	IF (ALLOCATED(delsurfstor)) DEALLOCATE(delsurfstor)
2509	!
2510	END SUBROUTINE routing_simple_clear
2511	!
2512
2513	!! ================================================================================================================================
2514	!! SUBROUTINE : routing_simple_lake
2515	!!
2516	!>\BRIEF : This subroutine stores water in lakes so that it does not cycle through the runoff.
2517	!! For the moment it only works for endoheric lakes but I can be extended in the future.
2518	!!
2519	!! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir
2520	!! is based on a maximum lake evaporation rate (maxevap_lake). \n
2521	!!
2522	!! RECENT CHANGE(S): None
2523	!!
2524	!! MAIN OUTPUT VARIABLE(S):
2525	!!
2526	!! REFERENCES : None
2527	!!
2528	!! FLOWCHART :None
2529	!! \n
2530	!_ ================================================================================================================================
2531
2532	SUBROUTINE routing_simple_lake(nbpt, dt_routing, lakeinflow, humrel, contfrac, return_lakes)
2533
2534	USE grid, ONLY : area
2535
2536	IMPLICIT NONE
2537	!! 0 Variable and parameter description
2538	!! 0.1 Input variables
2539
2540	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
2541	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
2542	REAL(r_std), INTENT(in) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt)
2543	REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless)
2544	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1)
2545
2546	!! 0.2 Output variables
2547	REAL(r_std), INTENT(out) :: return_lakes(nbpt) !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
2548
2549	!! 0.3 Local variables
2550	INTEGER(i_std) :: ig !! Indices (unitless)
2551	REAL(r_std) :: refill !!
2552	REAL(r_std) :: total_area !! Sum of all the surfaces of the basins (m^2)
2553
2554	!_ ================================================================================================================================
2555
2556
2557	DO ig=1,nbpt
2558	!
2559	total_area = area(ig)*contfrac(ig)
2560	!
2561	lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig)
2562	!uptake in Kg/dt
2563	refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area)
2564	return_lakes(ig) = MIN(refill, lake_reservoir(ig))
2565	lake_reservoir(ig) = lake_reservoir(ig) - return_lakes(ig)
2566	!Return in Kg/m^2/dt
2567	IF ( doswamps ) THEN
2568	return_lakes(ig) = return_lakes(ig)/total_area
2569	ELSE
2570	return_lakes(ig) = zero
2571	ENDIF
2572	!
2573	! This is the volume of the lake scaled to the entire grid.
2574	! It would be batter to scale it to the size of the lake
2575	! but this information is not yet available.
2576	lake_diag(ig) = lake_reservoir(ig)/total_area
2577	!
2578	ENDDO
2579
2580	END SUBROUTINE routing_simple_lake
2581
2582	END MODULE routing_simple

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

Download in other formats: