Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90 @ 7252

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

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

Download in other formats: