Context Navigation

routing_simple.f90 @ 8548

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

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

Context Navigation

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

Download in other formats: