Context Navigation

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

Last change on this file since 7485 was 7481, checked in by josefine.ghattas, 2 years ago
Enhencement on routing simple. Done and tested by Yann Meurdesoif and Arnaud Caubel
Property svn:keywords set to `Date Revision HeadURL`
File size: 139.2 KB

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

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

Download in other formats: