Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/routing.f90 @ 366

Last change on this file since 366 was 354, checked in by didier.solyga, 13 years ago
Add labels for the new externalized parameters. Replace 1000 by mille in the code
File size: 208.0 KB

Line
1	!!
2	!!
3	!! This module routes the water over the continents into the oceans.
4	!!
5	!! Histoire Salee
6	!!---------------
7	!! La douce riviere
8	!! Sortant de son lit
9	!! S'est jetee ma chere
10	!! dans les bras mais oui
11	!! du beau fleuve
12	!!
13	!! L'eau coule sous les ponts
14	!! Et puis les flots s'emeuvent
15	!! - N'etes vous pas au courant ?
16	!! Il parait que la riviere
17	!! Va devenir mer
18	!! Roland Bacri
19	!!
20	!! @author Jan Polcher
21	!! @Version : $Revision: 274 $, $Date: 2011-06-21 15:18:18 +0200 (Tue, 21 Jun 2011) $
22	!!
23	!< $HeadURL: http://forge.ipsl.jussieu.fr/orchidee/svn/trunk/ORCHIDEE/src_sechiba/routing.f90 $
24	!< $Date: 2011-06-21 15:18:18 +0200 (Tue, 21 Jun 2011) $
25	!< $Author: martial.mancip $
26	!< $Revision: 274 $
27	!! IPSL (2006)
28	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
29	!!
30	MODULE routing
31	!
32	!
33	! routines called : restput, restget
34	!
35	USE ioipsl
36	!
37	USE constantes
38	USE sechiba_io
39	USE pft_parameters
40	USE grid
41	USE parallel
42	!
43	IMPLICIT NONE
44	!
45	! public routines :
46	!
47	PRIVATE
48	PUBLIC :: routing_main, routing_clear
49	!
50	! variables used inside hydrol module : declaration and initialisation
51	!
52	LOGICAL, SAVE :: l_first_routing=.TRUE. !! Initialisation has to be done one time
53	LOGICAL, SAVE :: check_waterbal=.FALSE. !! The check the water balance
54	!
55	! The maximum number of basins we wish to have per grid-box (truncation of the model)
56	INTEGER(i_std), PARAMETER :: nbasmax=5
57	! The maximum number of bassins we can handle at any time during the generation of the maps.
58	INTEGER(i_std), PARAMETER :: nbvmax = 220
59	!
60	! The time constants are in days.
61	!
62	REAL(r_std), PARAMETER :: fast_tcst = 3.0, slow_tcst = 25.0, stream_tcst = 0.24
63	REAL(r_std), PARAMETER :: evap_cst = 0.18, wdelay_cst = 0.7
64	!
65	! Maximum evaporation rate from lakes 7.5 kg/m^2/d transformed into kg/m^2/sec
66	!
67	REAL(r_std), PARAMETER :: maxevap_lake = 7.5/86400.
68	!
69	REAL(r_std),SAVE :: dt_routing
70	!
71	! The variables describing the basins and their routing, need to be in the restart file.
72	!
73	INTEGER(i_std), SAVE :: num_largest
74	REAL(r_std), SAVE :: time_counter
75	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_loc
76	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_loc
77	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_loc
78	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_loc
79	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_loc
80	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_loc !! Variable to diagnose the hydrographs
81	!
82	! parallelism
83	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_glo
84	REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_glo
85	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_glo
86	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_glo
87	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_glo
88	INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_glo !! Variable to diagnose the hydrographs
89
90	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: routing_area
91	REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: topo_resid
92	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_togrid
93	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_tobasin
94	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: global_basinid
95	INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: hydrodiag
96	! Map of irrigated areas and floodplains
97	!
98	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigated ! Surface which can be irrigate in each grid-box
99	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodplains ! Surface which can be inondated in each grid-box
100	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: previous_outflow
101	!
102	! The reservoirs, also to be put into the restart file.
103	!
104	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_reservoir, slow_reservoir, stream_reservoir
105	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_reservoir
106	!
107	! The accumulated fluxes.
108	!
109	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: runoff_mean, drainage_mean, evapot_mean
110	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: precip_mean, humrel_mean, totnobio_mean, vegtot_mean
111	!
112	! The averaged outflow fluxes.
113	!
114	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lakeinflow_mean
115	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: returnflow_mean, irrigation_mean
116	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean, coastalflow_mean
117	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodtemp
118	INTEGER(i_std), SAVE :: floodtemp_lev
119	!
120	! Diagnostic variables ... well sort of !
121	!
122	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_netereq
123	!
124	! Hydrographs at the outflow of the grid box for major basins
125	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrographs
126	! Diagnostics for the various reservoirs we use (Kg/m^2)
127	REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: fast_diag, slow_diag, stream_diag, lake_diag
128	!
129	CONTAINS
130	!
131	!---------------------------------------------------------------------------------
132	!
133	SUBROUTINE routing_main(kjit, nbpt, dtradia, ldrestart_read, ldrestart_write, index, &
134	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
135	& drainage, evapot, precip_rain, humrel, &
136	& stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
137	!
138	IMPLICIT NONE
139	!
140	! This module will route the runoff and drainage produced by the hydrol module. These two
141	! fluxes are provided in kg/m^2dt. The result of the routing are 3 fluxes :
142	! - riverflow : The water which flows out from the major rivers. The flux will be located
143	! on the continental grid but this should be a coastal point.
144	! - coastalflow : This is the water which flows in a disperse way into the ocean. Essentially these
145	! are the outflows from all the small rivers.
146	! - returnflow : This is the water which flows into a land-point. Typically rivers which end in
147	! the desert. This water will go back into the hydrol module to allow re-evaporation.
148	! - irrigation : This is water taken from the river reservoir and beeing put into the upper
149	! layers of the soil.
150	! The two first fluxes are in kg/dt and the last one is on kg/m^2dt.
151	!
152	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
153	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size
154	INTEGER(i_std),INTENT(in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
155	INTEGER(i_std),INTENT(in) :: hist2_id !! _history_ file 2 identifier
156	REAL(r_std), INTENT(in) :: dtradia !! Time step in seconds
157	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
158	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
159	INTEGER(i_std), INTENT(in) :: index(nbpt) ! Indeces of the points on the map
160	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
161	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
162	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
163	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box
164	REAL(r_std), INTENT(in) :: totfrac_nobio(nbpt) ! Total fraction of continental ice+lakes+...
165	REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) ! Maximum vegetation fraction. We want to have the
166	! part of the grid which can have some vegetation.
167	REAL(r_std), INTENT(in) :: runoff(nbpt) ! grid-point runoff
168	REAL(r_std), INTENT(in) :: drainage(nbpt) ! grid-point drainage
169	REAL(r_std), INTENT(in) :: evapot(nbpt) ! Potential evaporation
170	REAL(r_std), INTENT(in) :: precip_rain(nbpt) ! Rainfall needed for the irrigation formula
171	REAL(r_std), INTENT(in) :: humrel(nbpt,nvm) ! Soil moisture stress
172	REAL(r_std), INTENT(in) :: stempdiag(nbpt,nbdl)! Temperature profile in soil
173	!
174	REAL(r_std), INTENT(out) :: returnflow(nbpt) ! The water flow which returns to the grid box (kg/m^2 per dt)
175	REAL(r_std), INTENT(out) :: irrigation(nbpt) ! Irrigation flux (kg/m^2 per dt)
176	REAL(r_std), INTENT(out) :: riverflow(nbpt) ! Outflow of the major rivers
177	REAL(r_std), INTENT(out) :: coastalflow(nbpt) ! Outflow on coastal points by small basins
178	!
179	! LOCAL
180	!
181	CHARACTER(LEN=30) :: var_name
182	REAL(r_std), DIMENSION(1) :: tmp_day
183	REAL(r_std), DIMENSION(nbpt) :: return_lakes
184	INTEGER(i_std) :: ig, jv
185	LOGICAL, SAVE :: init_irrig=.FALSE., init_flood=.FALSE.
186	!
187	! do initialisation
188	!
189	IF (l_first_routing) THEN
190	!
191	! Here we will allocate the memory and get the fixed fields from the restart file.
192	! If the info is not found then we will compute the routing map.
193	!
194	CALL routing_init (kjit, nbpt, index, dtradia, returnflow, irrigation, &
195	& riverflow, coastalflow, stempdiag, rest_id)
196	routing_area => routing_area_loc
197	topo_resid => topo_resid_loc
198	route_togrid => route_togrid_loc
199	route_tobasin => route_tobasin_loc
200	global_basinid => global_basinid_loc
201	hydrodiag => hydrodiag_loc
202	!
203	! This routine computes the routing map if needed.
204	!
205	IF ( COUNT(route_togrid_glo .GE. undef_int) .GT. 0 ) THEN
206	CALL routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac)
207	ENDIF
208	!
209	! Do we have what we need if we want to do irrigation
210	!
211	IF ( doirrigation .OR. dofloodplains ) THEN
212
213	IF ( doirrigation ) THEN
214	IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) THEN
215	init_irrig = .TRUE.
216	ENDIF
217	ENDIF
218
219	IF ( dofloodplains ) THEN
220	IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) THEN
221	init_flood = .TRUE.
222	ENDIF
223	ENDIF
224
225	IF ( init_irrig .OR. init_flood ) THEN
226	CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, &
227	contfrac, init_irrig, irrigated, init_flood, floodplains, hist_id, hist2_id)
228	ENDIF
229
230	ENDIF
231	!
232	! This routine gives a diag nostic of the basins used.
233	!
234	CALL routing_diagnostic_p(nbpt, index, resolution, contfrac, hist_id, hist2_id)
235	!
236	l_first_routing = .FALSE.
237	!
238	RETURN
239	!
240	ENDIF
241	!
242	! Accumulate the incoming water fluxes
243	!
244	runoff_mean(:) = runoff_mean(:) + runoff(:)
245	drainage_mean(:) = drainage_mean(:) + drainage(:)
246	evapot_mean(:) = evapot_mean(:) + evapot(:)
247	floodtemp(:) = stempdiag(:,floodtemp_lev)
248	precip_mean(:) = precip_mean(:) + precip_rain(:)
249	!
250	! Averaged variables (i.e. *dtradia/dt_routing)
251	!
252	totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dtradia/dt_routing
253	!
254	! Only potentially vegetated surfaces are taken into account. At the start of
255	! the growing seasons (when veget and veget_max are still different) we will
256	! more weight to these areas.
257	!
258	DO jv=2,nvm
259	DO ig=1,nbpt
260	humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)veget_max(ig,jv)dtradia/dt_routing
261	vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dtradia/dt_routing
262	ENDDO
263	ENDDO
264	!
265	time_counter = time_counter + dtradia
266	!
267	! If the time has come we do the routing.
268	!
269	IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN
270	!
271	! Check the water balance if needed
272	!
273	IF ( check_waterbal ) THEN
274	CALL routing_waterbal(nbpt, .TRUE., runoff_mean, drainage_mean, returnflow_mean, &
275	& irrigation_mean, riverflow_mean, coastalflow_mean)
276	ENDIF
277	!
278	! Make sure we do not flood north of 49N
279	!
280	DO ig=1,nbpt
281	IF ( lalo(ig,1) > 49.0 ) THEN
282	floodtemp(ig) = tp_00 - un
283	ENDIF
284	ENDDO
285	!
286	CALL routing_flow(nbpt, dt_routing, runoff_mean, drainage_mean, &
287	& vegtot_mean, totnobio_mean, evapot_mean, precip_mean, humrel_mean, floodtemp, &
288	& lakeinflow_mean, returnflow_mean, irrigation_mean, riverflow_mean, &
289	& coastalflow_mean, hydrographs)
290	!
291	CALL routing_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, return_lakes)
292	!
293	returnflow_mean(:) = returnflow_mean(:) + return_lakes(:)
294	!
295	! Close the water balance if asked for
296	!
297	IF ( check_waterbal ) THEN
298	CALL routing_waterbal(nbpt, .FALSE., runoff_mean, drainage_mean, returnflow_mean, &
299	& irrigation_mean, riverflow_mean, coastalflow_mean)
300	ENDIF
301	!
302	time_counter = zero
303	!
304	runoff_mean(:) = zero
305	drainage_mean(:) = zero
306	evapot_mean(:) = zero
307	precip_mean(:) = zero
308	!
309	humrel_mean(:) = zero
310	totnobio_mean(:) = zero
311	vegtot_mean(:) = zero
312	!
313	! Change the units of the routing fluxes from kg/dt_routing into kg/dtradia
314	! and from m^3/dt_routing into m^3/dtradia
315	!
316	returnflow_mean(:) = returnflow_mean(:)/dt_routing*dtradia
317	irrigation_mean(:) = irrigation_mean(:)/dt_routing*dtradia
318	irrig_netereq(:) = irrig_netereq(:)/dt_routing*dtradia
319	!
320	!
321	riverflow_mean(:) = riverflow_mean(:)/dt_routing*dtradia
322	coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dtradia
323	hydrographs(:) = hydrographs(:)/dt_routing*dtradia
324	!
325	! Convert from kg/dtradia to m^3/dtradia
326	!
327	hydrographs(:) = hydrographs(:)/mille
328	!
329	ENDIF
330	!
331	! Return the fraction of routed water for this time step.
332	!
333	returnflow(:) = returnflow_mean(:)
334	irrigation(:) = irrigation_mean(:)
335	riverflow(:) = riverflow_mean(:)
336	coastalflow(:) = coastalflow_mean(:)
337	!
338	! Write restart
339	!
340	IF (ldrestart_write) THEN
341	!
342	var_name ="routingcounter"
343	tmp_day(1) = time_counter
344	IF (is_root_prc) CALL restput (rest_id, var_name, 1, 1, 1, kjit, tmp_day)
345	!
346	var_name = 'routingarea'
347	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, routing_area, 'scatter', nbp_glo, index_g)
348	var_name = 'routetogrid'
349	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, REAL(route_togrid,r_std), 'scatter', &
350	& nbp_glo, index_g)
351	var_name = 'routetobasin'
352	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, REAL(route_tobasin,r_std), 'scatter', &
353	& nbp_glo, index_g)
354	var_name = 'basinid'
355	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, REAL(global_basinid,r_std), 'scatter', &
356	& nbp_glo, index_g)
357	var_name = 'topoindex'
358	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, topo_resid, 'scatter', nbp_glo, index_g)
359	var_name = 'fastres'
360	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, fast_reservoir, 'scatter', nbp_glo, index_g)
361	var_name = 'slowres'
362	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, slow_reservoir, 'scatter', nbp_glo, index_g)
363	var_name = 'streamres'
364	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, stream_reservoir, 'scatter',nbp_glo,index_g)
365	!
366	var_name = 'lakeres'
367	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, lake_reservoir, 'scatter', nbp_glo, index_g)
368	!
369	var_name = 'lakeinflow'
370	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g)
371	var_name = 'returnflow'
372	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, returnflow_mean, 'scatter', nbp_glo, index_g)
373	var_name = 'riverflow'
374	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g)
375	var_name = 'coastalflow'
376	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g)
377	var_name = 'hydrographs'
378	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, hydrographs, 'scatter', nbp_glo, index_g)
379	!
380	! Keep track of the accumulated variables
381	!
382	var_name = 'runoff_route'
383	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g)
384	var_name = 'drainage_route'
385	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g)
386	var_name = 'evapot_route'
387	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, evapot_mean, 'scatter', nbp_glo, index_g)
388	var_name = 'precip_route'
389	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, precip_mean, 'scatter', nbp_glo, index_g)
390	var_name = 'humrel_route'
391	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, humrel_mean, 'scatter', nbp_glo, index_g)
392	var_name = 'totnobio_route'
393	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, totnobio_mean, 'scatter', nbp_glo, index_g)
394	var_name = 'vegtot_route'
395	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, vegtot_mean, 'scatter', nbp_glo, index_g)
396	!
397	var_name = 'previous_outflow'
398	CALL restput_p (rest_id, var_name, nbp_glo, nbasmax+3, 1, kjit,previous_outflow,'scatter',nbp_glo,index_g)
399	!
400	IF ( doirrigation ) THEN
401	var_name = 'irrigated'
402	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, irrigated, 'scatter', nbp_glo, index_g)
403	var_name = 'irrigation'
404	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter', nbp_glo, index_g)
405	ENDIF
406	!
407	IF ( dofloodplains ) THEN
408	var_name = 'floodplains'
409	CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, floodplains, 'scatter', nbp_glo, index_g)
410	ENDIF
411	!
412	RETURN
413	!
414	ENDIF
415	!
416	! Write diagnostics
417	!
418	IF ( .NOT. almaoutput ) THEN
419	!
420	CALL histwrite(hist_id, 'riversret', kjit, returnflow, nbpt, index)
421	CALL histwrite(hist_id, 'hydrographs', kjit, hydrographs, nbpt, index)
422	!
423	CALL histwrite(hist_id, 'fastr', kjit, fast_diag, nbpt, index)
424	CALL histwrite(hist_id, 'slowr', kjit, slow_diag, nbpt, index)
425	CALL histwrite(hist_id, 'streamr', kjit, stream_diag, nbpt, index)
426	CALL histwrite(hist_id, 'lakevol', kjit, lake_diag, nbpt, index)
427	!
428	CALL histwrite(hist_id, 'irrigation', kjit, irrigation, nbpt, index)
429	!
430	CALL histwrite(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
431	!
432	ELSE
433	!
434	CALL histwrite(hist_id, 'dis', kjit, hydrographs, nbpt, index)
435	!
436	ENDIF
437	IF ( hist2_id > 0 ) THEN
438	IF ( .NOT. almaoutput ) THEN
439	!
440	CALL histwrite(hist2_id, 'riversret', kjit, returnflow, nbpt, index)
441	CALL histwrite(hist2_id, 'hydrographs', kjit, hydrographs, nbpt, index)
442	!
443	CALL histwrite(hist2_id, 'fastr', kjit, fast_diag, nbpt, index)
444	CALL histwrite(hist2_id, 'slowr', kjit, slow_diag, nbpt, index)
445	CALL histwrite(hist2_id, 'streamr', kjit, stream_diag, nbpt, index)
446	CALL histwrite(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index)
447	!
448	CALL histwrite(hist2_id, 'irrigation', kjit, irrigation, nbpt, index)
449	!
450	CALL histwrite(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
451	!
452	ELSE
453	!
454	CALL histwrite(hist2_id, 'dis', kjit, hydrographs, nbpt, index)
455	!
456	ENDIF
457	ENDIF
458	!
459	!
460	END SUBROUTINE routing_main
461	!
462	!---------------------------------------------------------------------------------
463	!
464	SUBROUTINE routing_init(kjit, nbpt, index, dtradia, returnflow, irrigation, &
465	& riverflow, coastalflow, stempdiag, rest_id)
466	!
467	IMPLICIT NONE
468	!
469	! interface description
470	! input
471	!
472	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
473	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size
474	INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index !! Indeces of the points on the map
475	REAL(r_std), INTENT(in) :: dtradia !! timestep
476	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: returnflow !! The water flow which returns to the grid box (kg/m^2 per dt)
477	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: irrigation !! Irrigation flow (kg/m^2 per dt)
478	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: riverflow !! Outflow of the major rivers
479	REAL(r_std), DIMENSION (nbpt),INTENT(out) :: coastalflow !! Outflow on coastal points by small basins
480	REAL(r_std), DIMENSION(nbpt,nbdl),INTENT(in) :: stempdiag !! Temperature profile in soil
481	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier
482	!
483	! local declaration
484	!
485	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
486	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g !! A temporary real array for the integers
487	REAL(r_std), DIMENSION(1) :: tmp_day
488	REAL(r_std) :: ratio, totarea
489	INTEGER(i_std) :: ier, ig, ib, ipn(1)
490	!
491	! These variables will require the configuration infrastructure
492	!
493	!Config Key = ROUTING_TIMESTEP
494	!Config If = RIVER_ROUTING
495	!Config Desc = Time step of th routing scheme
496	!Config Def = one_day
497	!Config Help = This values gives the time step in seconds of the routing scheme.
498	!Config It should be multiple of the main time step of ORCHIDEE. One day
499	!Config is a good value.
500	!
501	dt_routing = one_day
502	CALL getin_p('ROUTING_TIMESTEP', dt_routing)
503	!
504	!Config Key = ROUTING_RIVERS
505	!Config If = RIVER_ROUTING
506	!Config Desc = Number of rivers
507	!Config Def = 50
508	!Config Help = This parameter chooses the number of largest river basins
509	!Config which should be treated as independently as rivers and not
510	!Config flow into the oceans as diffusion coastal flow.
511	num_largest = 50
512	CALL getin_p('ROUTING_RIVERS', num_largest)
513	!
514	! In order to simplify the time cascade check that dt_routing
515	! is a multiple of dtradia
516	!
517	ratio = dt_routing/dtradia
518	IF ( ABS(NINT(ratio) - ratio) .GT. 10*EPSILON(ratio)) THEN
519	WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
520	WRITE(numout,*) "The chosen time step for the routing is not a multiple of the"
521	WRITE(numout,*) "main time step of the model. We will change dt_routing so that"
522	WRITE(numout,*) "this condition os fulfilled"
523	dt_routing = NINT(ratio) * dtradia
524	WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
525	ENDIF
526	!
527	IF ( dt_routing .LT. dtradia) THEN
528	WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
529	WRITE(numout,*) 'The routing timestep can not be smaller than the one'
530	WRITE(numout,*) 'of the model. We reset its value to the model''s timestep.'
531	dt_routing = dtradia
532	WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
533	ENDIF
534	!
535	var_name ="routingcounter"
536	IF (is_root_prc) THEN
537	CALL ioconf_setatt('UNITS', 's')
538	CALL ioconf_setatt('LONG_NAME','Time counter for the routing scheme')
539	CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day)
540	time_counter = tmp_day(1)
541	CALL setvar (time_counter, val_exp, 'NO_KEYWORD', zero)
542	ENDIF
543	CALL bcast(time_counter)
544	!!$ CALL setvar_p (time_counter, val_exp, 'NO_KEYWORD', zero)
545
546	!
547	ALLOCATE (routing_area_loc(nbpt,nbasmax))
548	ALLOCATE (routing_area_glo(nbp_glo,nbasmax))
549	var_name = 'routingarea'
550	IF (is_root_prc) THEN
551	CALL ioconf_setatt('UNITS', 'm^2')
552	CALL ioconf_setatt('LONG_NAME','Area of basin')
553	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g)
554	ENDIF
555	CALL scatter(routing_area_glo,routing_area_loc)
556	routing_area=>routing_area_loc
557	!
558	ALLOCATE (tmp_real_g(nbp_glo,nbasmax))
559	!
560	ALLOCATE (route_togrid_loc(nbpt,nbasmax))
561	ALLOCATE (route_togrid_glo(nbp_glo,nbasmax)) ! used in global in routing_flow
562	IF (is_root_prc) THEN
563	var_name = 'routetogrid'
564	CALL ioconf_setatt('UNITS', '-')
565	CALL ioconf_setatt('LONG_NAME','Grid into which the basin flows')
566	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
567	route_togrid_glo(:,:) = undef_int
568	WHERE ( tmp_real_g .LT. val_exp )
569	route_togrid_glo = NINT(tmp_real_g)
570	ENDWHERE
571	ENDIF
572	CALL bcast(route_togrid_glo) ! used in global in routing_flow
573	CALL scatter(route_togrid_glo,route_togrid_loc)
574	route_togrid=>route_togrid_loc
575	!
576	ALLOCATE (route_tobasin_loc(nbpt,nbasmax))
577	ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax))
578	IF (is_root_prc) THEN
579	var_name = 'routetobasin'
580	CALL ioconf_setatt('UNITS', '-')
581	CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes')
582	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
583	route_tobasin_glo = undef_int
584	WHERE ( tmp_real_g .LT. val_exp )
585	route_tobasin_glo = NINT(tmp_real_g)
586	ENDWHERE
587	ENDIF
588	CALL scatter(route_tobasin_glo,route_tobasin_loc)
589	route_tobasin=>route_tobasin_loc
590	!
591	ALLOCATE (global_basinid_loc(nbpt,nbasmax))
592	ALLOCATE (global_basinid_glo(nbp_glo,nbasmax))
593	IF (is_root_prc) THEN
594	var_name = 'basinid'
595	CALL ioconf_setatt('UNITS', '-')
596	CALL ioconf_setatt('LONG_NAME','ID of basin')
597	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
598	global_basinid_glo = undef_int
599	WHERE ( tmp_real_g .LT. val_exp )
600	global_basinid_glo = NINT(tmp_real_g)
601	ENDWHERE
602	ENDIF
603	CALL scatter(global_basinid_glo,global_basinid_loc)
604	global_basinid=>global_basinid_loc
605	!
606	ALLOCATE (topo_resid_loc(nbpt,nbasmax))
607	ALLOCATE (topo_resid_glo(nbp_glo,nbasmax))
608	IF (is_root_prc) THEN
609	var_name = 'topoindex'
610	CALL ioconf_setatt('UNITS', 'm')
611	CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time')
612	CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g)
613	ENDIF
614	CALL scatter(topo_resid_glo,topo_resid_loc)
615	topo_resid=>topo_resid_loc
616	!
617	ALLOCATE (fast_reservoir(nbpt,nbasmax))
618	var_name = 'fastres'
619	CALL ioconf_setatt('UNITS', 'Kg')
620	CALL ioconf_setatt('LONG_NAME','Water in the fast reservoir')
621	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g)
622	CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero)
623	!
624	ALLOCATE (slow_reservoir(nbpt,nbasmax))
625	var_name = 'slowres'
626	CALL ioconf_setatt('UNITS', 'Kg')
627	CALL ioconf_setatt('LONG_NAME','Water in the slow reservoir')
628	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g)
629	CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero)
630	!
631	ALLOCATE (stream_reservoir(nbpt,nbasmax))
632	var_name = 'streamres'
633	CALL ioconf_setatt('UNITS', 'Kg')
634	CALL ioconf_setatt('LONG_NAME','Water in the stream reservoir')
635	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g)
636	CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero)
637	!
638	ALLOCATE (lake_reservoir(nbpt))
639	var_name = 'lakeres'
640	CALL ioconf_setatt('UNITS', 'Kg')
641	CALL ioconf_setatt('LONG_NAME','Water in the lake reservoir')
642	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g)
643	CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero)
644	!
645	! Map of irrigated areas
646	!
647	IF ( doirrigation ) THEN
648	ALLOCATE (irrigated(nbpt))
649	var_name = 'irrigated'
650	CALL ioconf_setatt('UNITS', 'm^2')
651	CALL ioconf_setatt('LONG_NAME','Surface of irrigated area')
652	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g)
653	CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba)
654	ENDIF
655	!
656	ALLOCATE (previous_outflow(nbpt, nbasmax+3))
657	var_name = 'previous_outflow'
658	CALL ioconf_setatt('UNITS', 'm^3/dt')
659	CALL ioconf_setatt('LONG_NAME','Previous outflow from this basin')
660	CALL restget_p (rest_id, var_name, nbp_glo, nbasmax+3, 1, kjit, .TRUE., previous_outflow, "gather", nbp_glo, index_g)
661	CALL setvar_p (previous_outflow, val_exp, 'NO_KEYWORD', zero)
662	!
663	IF ( dofloodplains ) THEN
664	ALLOCATE (floodplains(nbpt))
665	var_name = 'floodplains'
666	CALL ioconf_setatt('UNITS', 'm^2')
667	CALL ioconf_setatt('LONG_NAME','Surface which can be flooded')
668	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodplains, "gather", nbp_glo, index_g)
669	CALL setvar_p (floodplains, val_exp, 'NO_KEYWORD', undef_sechiba)
670	ENDIF
671	!
672	! Put into the restart file the fluxes so that they can be regenerated at restart.
673	!
674	ALLOCATE (lakeinflow_mean(nbpt))
675	var_name = 'lakeinflow'
676	CALL ioconf_setatt('UNITS', 'Kg/dt')
677	CALL ioconf_setatt('LONG_NAME','Lake inflow')
678	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
679	CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
680	!
681	ALLOCATE (returnflow_mean(nbpt))
682	var_name = 'returnflow'
683	CALL ioconf_setatt('UNITS', 'Kg/m^2/dt')
684	CALL ioconf_setatt('LONG_NAME','Deep return flux')
685	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g)
686	CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero)
687	returnflow(:) = returnflow_mean(:)
688	!
689	ALLOCATE (irrigation_mean(nbpt))
690	ALLOCATE (irrig_netereq(nbpt))
691	irrig_netereq(:) = zero
692	!
693	IF ( doirrigation ) THEN
694	var_name = 'irrigation'
695	CALL ioconf_setatt('UNITS', 'Kg/dt')
696	CALL ioconf_setatt('LONG_NAME','Artificial irrigation flux')
697	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g)
698	CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero)
699	ELSE
700	irrigation_mean(:) = zero
701	ENDIF
702	irrigation(:) = irrigation_mean(:)
703	!
704	ALLOCATE (riverflow_mean(nbpt))
705	var_name = 'riverflow'
706	CALL ioconf_setatt('UNITS', 'Kg/dt')
707	CALL ioconf_setatt('LONG_NAME','River flux into the sea')
708	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
709	CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
710	riverflow(:) = riverflow_mean(:)
711	!
712	ALLOCATE (coastalflow_mean(nbpt))
713	var_name = 'coastalflow'
714	CALL ioconf_setatt('UNITS', 'Kg/dt')
715	CALL ioconf_setatt('LONG_NAME','Diffuse flux into the sea')
716	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
717	CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
718	coastalflow(:) = coastalflow_mean(:)
719	!
720	! Locate it at the 2m level
721	ipn = MINLOC(ABS(diaglev-2))
722	floodtemp_lev = ipn(1)
723	ALLOCATE (floodtemp(nbpt))
724	floodtemp(:) = stempdiag(:,floodtemp_lev)
725	!
726	ALLOCATE(hydrographs(nbpt))
727	var_name = 'hydrographs'
728	CALL ioconf_setatt('UNITS', 'm^3/dt')
729	CALL ioconf_setatt('LONG_NAME','Hydrograph at outlow of grid')
730	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g)
731	CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero)
732	!
733	! The diagnostic variables, they are initialized from the above restart variables.
734	!
735	ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), lake_diag(nbpt), stat=ier)
736	!
737	fast_diag(:) = zero
738	slow_diag(:) = zero
739	stream_diag(:) = zero
740	lake_diag(:) = zero
741	!
742	DO ig=1,nbpt
743	totarea = zero
744	DO ib=1,nbasmax
745	totarea = totarea + routing_area(ig,ib)
746	fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib)
747	slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib)
748	stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib)
749	ENDDO
750	!
751	fast_diag(ig) = fast_diag(ig)/totarea
752	slow_diag(ig) = slow_diag(ig)/totarea
753	stream_diag(ig) = stream_diag(ig)/totarea
754	!
755	! This is the volume of the lake scaled to the entire grid.
756	! It would be batter to scale it to the size of the lake
757	! but this information is not yet available.
758	!
759	lake_diag(ig) = lake_reservoir(ig)/totarea
760	!
761	ENDDO
762	!
763	!
764	! Get from the restart the fluxes we accumulated.
765	!
766	ALLOCATE (runoff_mean(nbpt))
767	var_name = 'runoff_route'
768	CALL ioconf_setatt('UNITS', 'Kg')
769	CALL ioconf_setatt('LONG_NAME','Accumulated runoff for routing')
770	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
771	CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
772	!
773	ALLOCATE(drainage_mean(nbpt))
774	var_name = 'drainage_route'
775	CALL ioconf_setatt('UNITS', 'Kg')
776	CALL ioconf_setatt('LONG_NAME','Accumulated drainage for routing')
777	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
778	CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
779	!
780	ALLOCATE(evapot_mean(nbpt))
781	var_name = 'evapot_route'
782	CALL ioconf_setatt('UNITS', 'Kg/m^2')
783	CALL ioconf_setatt('LONG_NAME','Accumulated potential evaporation for routing')
784	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot_mean, "gather", nbp_glo, index_g)
785	CALL setvar_p (evapot_mean, val_exp, 'NO_KEYWORD', zero)
786	!
787	ALLOCATE(precip_mean(nbpt))
788	var_name = 'precip_route'
789	CALL ioconf_setatt('UNITS', 'Kg/m^2')
790	CALL ioconf_setatt('LONG_NAME','Accumulated rain precipitation for irrigation')
791	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g)
792	CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero)
793	!
794	ALLOCATE(humrel_mean(nbpt))
795	var_name = 'humrel_route'
796	CALL ioconf_setatt('UNITS', '-')
797	CALL ioconf_setatt('LONG_NAME','Mean humrel for irrigation')
798	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g)
799	CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un)
800	!
801	ALLOCATE(totnobio_mean(nbpt))
802	var_name = 'totnobio_route'
803	CALL ioconf_setatt('UNITS', '-')
804	CALL ioconf_setatt('LONG_NAME','Last Total fraction of no bio for irrigation')
805	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g)
806	CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero)
807	!
808	ALLOCATE(vegtot_mean(nbpt))
809	var_name = 'vegtot_route'
810	CALL ioconf_setatt('UNITS', '-')
811	CALL ioconf_setatt('LONG_NAME','Last Total fraction of vegetation')
812	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g)
813	CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un)
814	!
815	!
816	DEALLOCATE(tmp_real_g)
817	!
818	! Allocate diagnostic variables
819	!
820	ALLOCATE(hydrodiag_loc(nbpt,nbasmax),hydrodiag_glo(nbp_glo,nbasmax),stat=ier)
821	IF (ier.NE.0) THEN
822	WRITE (numout,) ' error in hydrodiag allocation. We stop. We need kjpindex words = ', nbptnbasmax
823	STOP 'routing_init'
824	END IF
825
826	hydrodiag=>hydrodiag_loc
827	!
828	!
829	END SUBROUTINE routing_init
830	!
831	!---------------------------------------------------------------------------------
832	!
833	SUBROUTINE routing_clear()
834	!
835	l_first_routing=.TRUE.
836	!
837	IF (ALLOCATED(routing_area_loc)) DEALLOCATE(routing_area_loc)
838	IF (ALLOCATED(route_togrid_loc)) DEALLOCATE(route_togrid_loc)
839	IF (ALLOCATED(route_tobasin_loc)) DEALLOCATE(route_tobasin_loc)
840	IF (ALLOCATED(global_basinid_loc)) DEALLOCATE(global_basinid_loc)
841	IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc)
842	IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo)
843	IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo)
844	IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo)
845	IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo)
846	IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo)
847	IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir)
848	IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir)
849	IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir)
850	IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir)
851	IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_mean)
852	IF (ALLOCATED(riverflow_mean)) DEALLOCATE(riverflow_mean)
853	IF (ALLOCATED(coastalflow_mean)) DEALLOCATE(coastalflow_mean)
854	IF (ALLOCATED(lakeinflow_mean)) DEALLOCATE(lakeinflow_mean)
855	IF (ALLOCATED(runoff_mean)) DEALLOCATE(runoff_mean)
856	IF (ALLOCATED(drainage_mean)) DEALLOCATE(drainage_mean)
857	IF (ALLOCATED(evapot_mean)) DEALLOCATE(evapot_mean)
858	IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean)
859	IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean)
860	IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean)
861	IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean)
862	IF (ALLOCATED(floodtemp)) DEALLOCATE(floodtemp)
863	IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc)
864	IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo)
865	IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs)
866	IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean)
867	IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated)
868	IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains)
869	!
870	END SUBROUTINE routing_clear
871	!
872	!---------------------------------------------------------------------------------
873	!
874	SUBROUTINE routing_flow(nbpt, dt_routing, runoff, drainage, &
875	& vegtot, totnobio, evapot, precip, humrel, floodtemp, &
876	& lakeinflow, returnflow, irrigation, riverflow, &
877	& coastalflow, hydrographs)
878	!
879	IMPLICIT NONE
880	!
881	!
882	! INPUT
883	!
884	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size
885	REAL(r_std), INTENT (in) :: dt_routing !! Time step in seconds
886	REAL(r_std), INTENT(in) :: runoff(nbpt) !! grid-point runoff
887	REAL(r_std), INTENT(in) :: drainage(nbpt) !! grid-point drainage
888	REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated area
889	REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas whichcan not have vegetation
890	REAL(r_std), INTENT(in) :: evapot(nbpt) !! grid-point potential evaporation
891	REAL(r_std), INTENT(in) :: precip(nbpt) !! Precipiation (rainfall here)
892	REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress
893	REAL(r_std), INTENT(in) :: floodtemp(nbpt) !! Temperature to decide if floodplains work
894	REAL(r_std), INTENT(out) :: lakeinflow(nbpt) !! The water flow which flows into lakes (kg/dt)
895	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! Water flowing back into soil moisture (kg/m^2/dt)
896	REAL(r_std), INTENT(out) :: irrigation(nbpt) !! The artificial irrigation (kg/m^2 per dt)
897	REAL(r_std), INTENT(out) :: riverflow(nbpt) !! Outflow of the major rivers (kg/dt)
898	REAL(r_std), INTENT(out) :: coastalflow(nbpt) !! Outflow on coastal points by small basins (kg/dt)
899	REAL(r_std), INTENT(out) :: hydrographs(nbpt) !! Hydrographs at the outflow of the gird box for major basins
900	!
901	! LOCAL
902	!
903	REAL(r_std), DIMENSION(nbpt, nbasmax) :: fast_flow
904	REAL(r_std), DIMENSION(nbpt, nbasmax) :: slow_flow
905	REAL(r_std), DIMENSION(nbpt, nbasmax) :: stream_flow
906	REAL(r_std), DIMENSION(nbpt, nbasmax) :: reinfiltration
907	REAL(r_std), DIMENSION(nbpt, nbasmax) :: baseirrig !! Irrigation uptake from each basin reservoir.
908	REAL(r_std), DIMENSION(nbpt, 0:nbasmax+3) :: transport
909	REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_glo
910	REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_sum
911	REAL(r_std), DIMENSION(nbpt, nbasmax) :: floods, wdelay, previous_outflow, potflood, inflow
912	REAL(r_std), DIMENSION(nbpt) :: tobeflooded
913	REAL(r_std), DIMENSION(nbpt) :: totarea
914	REAL(r_std) :: flow, floodindex
915	INTEGER(i_std) :: ig, ib, rtg, rtb, ierr
916
917	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_flow_g,slow_flow_g,stream_flow_g,floods_g,wdelay_g
918
919	!
920	transport(:,:) = zero
921	transport_glo(:,:) = zero
922	previous_outflow(:,:) = zero
923	irrig_netereq(:) = zero
924	baseirrig(:,:) = zero
925	totarea(:) = zero
926	!
927	! Compute all the fluxes
928	!
929	DO ig=1,nbpt
930	DO ib=1,nbasmax
931	!
932	totarea(ig) = totarea(ig) + routing_area(ig,ib)
933	!
934	IF ( route_tobasin(ig,ib) .GT. 0 ) THEN
935	!
936	flow = MIN(fast_reservoir(ig,ib)/((topo_resid(ig,ib)/mille)fast_tcstone_day/dt_routing),&
937	& fast_reservoir(ig,ib))
938	fast_flow(ig,ib) = flow
939	!
940	flow = MIN(slow_reservoir(ig,ib)/((topo_resid(ig,ib)/mille)slow_tcstone_day/dt_routing),&
941	& slow_reservoir(ig,ib))
942	slow_flow(ig,ib) = flow
943	!
944	flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/mille)stream_tcstone_day/dt_routing),&
945	& stream_reservoir(ig,ib))
946	stream_flow(ig,ib) = flow
947	!
948	!
949	ELSE
950	fast_flow(ig,ib) = zero
951	slow_flow(ig,ib) = zero
952	stream_flow(ig,ib) = zero
953	ENDIF
954	inflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib)
955	ENDDO
956	ENDDO
957	!
958	! Do the floodings
959	!
960	IF ( dofloodplains ) THEN
961	DO ig=1,nbpt
962	tobeflooded(ig) = floodplains(ig)
963	DO ib=1,nbasmax
964	potflood(ig,ib) = inflow(ig,ib) - previous_outflow(ig,ib)
965	!
966	IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > zero .AND. floodtemp(ig) > tp_00 ) THEN
967	!
968	IF (routing_area(ig,ib) > tobeflooded(ig)) THEN
969	floodindex = tobeflooded(ig) / routing_area(ig,ib)
970	ELSE
971	floodindex = un
972	ENDIF
973	!
974	floods(ig,ib) = floodindex * evap_cst * potflood(ig,ib)
975	!
976	wdelay(ig,ib) = floodindex * wdelay_cst * potflood(ig,ib)
977	!
978	!
979	tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib)
980	!
981	ELSE
982	floods(ig,ib) = zero
983	wdelay(ig,ib) = zero
984	ENDIF
985	ENDDO
986	ENDDO
987	ELSE
988	DO ig=1,nbpt
989	DO ib=1,nbasmax
990	floods(ig,ib) = zero
991	wdelay(ig,ib) = zero
992	ENDDO
993	ENDDO
994	ENDIF
995	!
996	DO ig=1,nbpt
997	DO ib=1,nbasmax
998	reinfiltration(ig,ib) = floods(ig,ib)
999	ENDDO
1000	ENDDO
1001
1002	!ym cette methode conserve les erreurs d'arrondie
1003	!ym mais n'est pas la plus efficace
1004
1005	IF (is_root_prc) THEN
1006	ALLOCATE( fast_flow_g(nbp_glo, nbasmax), slow_flow_g(nbp_glo, nbasmax), &
1007	stream_flow_g(nbp_glo, nbasmax), floods_g(nbp_glo, nbasmax), &
1008	wdelay_g(nbp_glo, nbasmax) )
1009	ELSE
1010	ALLOCATE( fast_flow_g(1,1), slow_flow_g(1,1), &
1011	stream_flow_g(1, 1), floods_g(1,1), &
1012	wdelay_g(1,1) )
1013	ENDIF
1014
1015
1016	CALL gather(fast_flow,fast_flow_g)
1017	CALL gather(slow_flow,slow_flow_g)
1018	CALL gather(stream_flow,stream_flow_g)
1019	CALL gather(floods,floods_g)
1020	CALL gather(wdelay,wdelay_g)
1021
1022	IF (is_root_prc) THEN
1023	DO ig=1,nbp_glo
1024	DO ib=1,nbasmax
1025	!
1026	rtg = route_togrid_glo(ig,ib)
1027	rtb = route_tobasin_glo(ig,ib)
1028	transport_glo(rtg,rtb) = transport_glo(rtg,rtb) + fast_flow_g(ig,ib) + slow_flow_g(ig,ib) + &
1029	& stream_flow_g(ig,ib) - floods_g(ig,ib) - wdelay_g(ig,ib)
1030	!
1031	ENDDO
1032	ENDDO
1033	ENDIF
1034
1035	DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g, floods_g, wdelay_g )
1036
1037	CALL scatter(transport_glo,transport)
1038
1039	!ym DO ig=1,nbpt
1040	!ym DO ib=1,nbasmax
1041	!ym !
1042	!ym rtg = route_togrid(ig,ib)
1043	!ym rtb = route_tobasin(ig,ib)
1044	!ym transport_glo(rtg,rtb) = transport_glo(rtg,rtb) + fast_flow(ig,ib) + slow_flow(ig,ib) + &
1045	!ym & stream_flow(ig,ib) - floods(ig,ib) - wdelay(ig,ib)
1046	!ym !
1047	!ym ENDDO
1048	!ym ENDDO
1049	!ym
1050	!ym CALL reduce_sum(transport_glo,transport_sum)
1051	!ym CALL scatter(transport_sum,transport)
1052
1053	!
1054	! Update all reservoirs
1055	!
1056	DO ig=1,nbpt
1057	DO ib=1,nbasmax
1058	fast_reservoir(ig,ib) = fast_reservoir(ig,ib) + runoff(ig)*routing_area(ig,ib) - &
1059	& reinfiltration(ig,ib) - fast_flow(ig,ib) + floods(ig,ib) + wdelay(ig,ib)
1060	slow_reservoir(ig,ib) = slow_reservoir(ig,ib) + drainage(ig)*routing_area(ig,ib) - &
1061	& slow_flow(ig,ib)
1062	stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + transport(ig,ib) - &
1063	& stream_flow(ig,ib)
1064	!
1065	previous_outflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) - &
1066	& wdelay(ig,ib) - floods(ig,ib)
1067	!
1068	ENDDO
1069	ENDDO
1070	!
1071	! Compute the total reinfiltration in the grid box
1072	!
1073	returnflow(:) = zero
1074	DO ig=1,nbpt
1075	!
1076	DO ib=1,nbasmax
1077	returnflow(ig) = returnflow(ig) + reinfiltration(ig,ib)
1078	ENDDO
1079	!
1080	returnflow(ig) = returnflow(ig)/totarea(ig)
1081	!
1082	ENDDO
1083	!
1084	! Compute the net irrigation requirement from Univ of Kassel
1085	!
1086	! This is a very low priority process and thus only applies if
1087	! there is some water left in the reservoirs after all other things.
1088	!
1089	IF ( doirrigation ) THEN
1090	DO ig=1,nbpt
1091
1092	IF ((vegtot(ig) .GT. zero) .AND. (humrel(ig) .LT. 0.99)) THEN
1093	irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, &
1094	& crop_coef * evapot(ig) - &
1095	& MAX(precip(ig)+returnflow(ig)-runoff(ig)-drainage(ig), zero) )
1096	irrig_netereq(ig) = 1 * irrig_netereq(ig)
1097
1098	IF(irrig_netereq(ig).LT.zero) THEN
1099	WRITE(numout,*) 'there is a probleme for irrig_netereq',ig,irrig_netereq(ig)
1100	ENDIF
1101
1102	ENDIF
1103
1104	DO ib=1,nbasmax
1105	IF ( route_tobasin(ig,ib) .GT. 0 ) THEN
1106
1107	baseirrig(ig,ib) = MIN( irrig_netereq(ig) * routing_area(ig,ib),&
1108	& stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) )
1109
1110	slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + &
1111	& MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-baseirrig(ig,ib))))
1112
1113	fast_reservoir(ig,ib) = MAX( zero, &
1114	& fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-baseirrig(ig,ib)))
1115	stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-baseirrig(ig,ib) )
1116
1117	IF(baseirrig(ig,ib) .LT. zero .OR. slow_reservoir(ig,ib) .LT. zero .OR. &
1118	& fast_reservoir(ig,ib) .LT. zero .OR. stream_reservoir(ig,ib) .LT. zero) THEN
1119	WRITE(numout,*) 'There is negative values related to irrigation', ig,ib,baseirrig(ig,ib), &
1120	& slow_reservoir(ig,ib),fast_reservoir(ig,ib),stream_reservoir(ig,ib)
1121	ENDIF
1122
1123	ENDIF
1124
1125	previous_outflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) - &
1126	& wdelay(ig,ib) - floods(ig,ib)
1127
1128	ENDDO
1129	ENDDO
1130	!
1131	ENDIF
1132	!
1133	!
1134	! Compute the fluxes which leave the routing scheme
1135	!
1136	! Lakeinflow is in Kg/dt
1137	! returnflow is in Kg/m^2/dt
1138	!
1139	hydrographs(:) = zero
1140	fast_diag(:) = zero
1141	slow_diag(:) = zero
1142	stream_diag(:) = zero
1143	irrigation(:) = zero
1144	!
1145	!
1146	DO ig=1,nbpt
1147	!
1148	DO ib=1,nbasmax
1149	IF (hydrodiag(ig,ib) > 0 ) THEN
1150	hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + &
1151	& stream_flow(ig,ib) - floods(ig,ib) - wdelay(ig,ib)
1152	ENDIF
1153	fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib)
1154	slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib)
1155	stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib)
1156	irrigation (ig) = irrigation (ig) + baseirrig(ig,ib)
1157	ENDDO
1158	!
1159	fast_diag(ig) = fast_diag(ig)/totarea(ig)
1160	slow_diag(ig) = slow_diag(ig)/totarea(ig)
1161	stream_diag(ig) = stream_diag(ig)/totarea(ig)
1162	!
1163	irrigation(ig) = irrigation(ig)/totarea(ig)
1164	!
1165	! The three output types for the routing : endoheric basins,, rivers and
1166	! diffuse coastal flow.
1167	!
1168	lakeinflow(ig) = transport(ig,nbasmax+1)
1169	coastalflow(ig) = transport(ig,nbasmax+2)
1170	riverflow(ig) = transport(ig,nbasmax+3)
1171	!
1172	IF ( irrigation(ig) .GE. irrig_netereq(ig)+1e4 ) THEN
1173	WRITE(numout,*) 'There is a problem here with irrigation',ig,irrigation(ig),irrig_netereq(ig)
1174	WRITE(numout,*) irrigated(ig),totarea(ig), evapot(ig), precip_mean(ig),runoff(ig),drainage(ig)
1175	STOP
1176	ENDIF
1177	!
1178	ENDDO
1179	!
1180	END SUBROUTINE routing_flow
1181	!
1182	!---------------------------------------------------------------------------------
1183	!
1184	SUBROUTINE routing_lake(nbpt, dt_routing, lakeinflow, humrel, returnflow)
1185	!
1186	IMPLICIT NONE
1187	!
1188	! This routine stores water in lakes so that it does not cycle through
1189	! the runoff. For the moment it only works for endoheric lakes but I can
1190	! be extended in the future.
1191	! The return flow to the soil moisture reservoir is based on a maximum
1192	! lake evaporation rate (maxevap_lake).
1193	!
1194	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size
1195	REAL(r_std), INTENT (in) :: dt_routing !! Time step in seconds
1196	REAL(r_std), INTENT(in) :: lakeinflow(nbpt) !! Flow into the lake (kg/dt)
1197	REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture deficit around the lake (Hum !)
1198	REAL(r_std), INTENT(out) :: returnflow(nbpt) !! Water flowing back into soil moisture (kg/m^2/dt)
1199	!
1200	! LOCAL
1201	!
1202	INTEGER(i_std) :: ig
1203	REAL(r_std) :: refill, total_area
1204	!
1205	!
1206	!
1207	DO ig=1,nbpt
1208	!
1209	total_area = SUM(routing_area(ig,:))
1210	!
1211	lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig)
1212	!uptake in Kg/dt
1213	refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area)
1214	returnflow(ig) = MIN(refill, lake_reservoir(ig))
1215	lake_reservoir(ig) = lake_reservoir(ig) - returnflow(ig)
1216	!Return in Kg/m^2/dt
1217	returnflow(ig) = returnflow(ig)/total_area
1218	!
1219	! This is the volume of the lake scaled to the entire grid.
1220	! It would be batter to scale it to the size of the lake
1221	! but this information is not yet available.
1222	lake_diag(ig) = lake_reservoir(ig)/total_area
1223	!
1224	ENDDO
1225	!
1226	END SUBROUTINE routing_lake
1227	!
1228	!---------------------------------------------------------------------------------
1229	!
1230	SUBROUTINE routing_diagnostic_p(nbpt,index, resolution, contfrac, hist_id, hist2_id)
1231	!
1232	IMPLICIT NONE
1233
1234	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size
1235	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indeces of the points on the map
1236	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
1237	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
1238	INTEGER(i_std),INTENT (in) :: hist_id !! _history_ file identifier
1239	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
1240	REAL(r_std), DIMENSION(nbpt) :: nbrivers !! Number of rivers in the grid
1241	REAL(r_std), DIMENSION(nbpt) :: basinmap !! Map of basins
1242	REAL(r_std), DIMENSION(nbp_glo) :: nbrivers_g !! Number of rivers in the grid
1243	REAL(r_std), DIMENSION(nbp_glo) :: basinmap_g !! Map of basins
1244
1245	routing_area => routing_area_glo
1246	topo_resid => topo_resid_glo
1247	route_togrid => route_togrid_glo
1248	route_tobasin => route_tobasin_glo
1249	global_basinid => global_basinid_glo
1250	hydrodiag=>hydrodiag_glo
1251
1252	IF (is_root_prc) CALL routing_diagnostic(nbp_glo,index_g, resolution_g, contfrac_g, nbrivers_g,basinmap_g)
1253
1254	routing_area => routing_area_loc
1255	topo_resid => topo_resid_loc
1256	route_togrid => route_togrid_loc
1257	route_tobasin => route_tobasin_loc
1258	global_basinid => global_basinid_loc
1259	hydrodiag=>hydrodiag_loc
1260
1261	CALL scatter(nbrivers_g,nbrivers)
1262	CALL scatter(basinmap_g,basinmap)
1263	CALL scatter(hydrodiag_glo,hydrodiag_loc)
1264
1265	IF ( .NOT. almaoutput ) THEN
1266	CALL histwrite(hist_id, 'basinmap', 1, basinmap, nbpt, index)
1267	CALL histwrite(hist_id, 'nbrivers', 1, nbrivers, nbpt, index)
1268	ELSE
1269	ENDIF
1270	IF ( hist2_id > 0 ) THEN
1271	IF ( .NOT. almaoutput ) THEN
1272	CALL histwrite(hist2_id, 'basinmap', 1, basinmap, nbpt, index)
1273	CALL histwrite(hist2_id, 'nbrivers', 1, nbrivers, nbpt, index)
1274	ELSE
1275	ENDIF
1276	ENDIF
1277
1278
1279	END SUBROUTINE routing_diagnostic_p
1280
1281
1282	SUBROUTINE routing_diagnostic(nbpt,index, resolution, contfrac,nbrivers,basinmap)
1283	!
1284	! This subroutine will set up a map of the major basins
1285	!
1286	! INPUTS
1287	!
1288	INTEGER(i_std), INTENT(in) :: nbpt !! Domain size
1289	INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indeces of the points on the map
1290	REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y
1291	REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box.
1292	!
1293	! OUTPUTS
1294	!
1295	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: nbrivers !! Number of rivers in the grid
1296	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: basinmap !! Map of basins
1297	!
1298	! LOCAL
1299	!
1300	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: pts !! list the points belonging to the basin
1301	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: ptbas !! list the basin number for this point
1302	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: outpt !! Outflow point for each basin
1303	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: nb_pts !! Number of points in the basin
1304	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: totarea !! Total area of basin
1305	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: topids !! The IDs of the first num_largest basins.
1306	CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: basin_names !! Names of the rivers
1307	CHARACTER(LEN=25) :: name_str
1308	!
1309	INTEGER(i_std) :: ig, ib, og, ob, ign, ibn, ff(1), ic, icc, nb_small
1310	!
1311	ALLOCATE(pts(num_largest, nbpt))
1312	ALLOCATE(ptbas(num_largest, nbpt))
1313	ALLOCATE(outpt(num_largest, 2))
1314	ALLOCATE(nb_pts(num_largest))
1315	ALLOCATE(totarea(num_largest))
1316	ALLOCATE(topids(num_largest))
1317	!
1318	!
1319	! First we get the list of all river outflow points
1320	! We work under the assumption that we only have num_largest basins finishing with
1321	! nbasmax+3. This is checked in routing_truncate.
1322	!
1323	nb_small = 1
1324	outpt(:,:) = -1
1325	ic = 0
1326	DO ig=1,nbpt
1327	DO ib=1,nbasmax
1328	ign = route_togrid(ig, ib)
1329	ibn = route_tobasin(ig, ib)
1330	IF ( ibn .EQ. nbasmax+3) THEN
1331	ic = ic + 1
1332	outpt(ic,1) = ig
1333	outpt(ic,2) = ib
1334	!
1335	! Get the largest id of the basins we call a river. This is
1336	! to extract the names of all rivers.
1337	!
1338	IF ( global_basinid(ig,ib) > nb_small ) THEN
1339	nb_small = global_basinid(ig,ib)
1340	ENDIF
1341	ENDIF
1342	ENDDO
1343	ENDDO
1344	!
1345	nb_small = MIN(nb_small, 349)
1346	!
1347	ALLOCATE(basin_names(nb_small))
1348	!
1349	CALL routing_names(nb_small, basin_names)
1350	!
1351	! Go through all points and basins to see if they outflow as a river and store the
1352	! information needed in the various arrays.
1353	!
1354	nb_pts(:) = 0
1355	totarea(:) = zero
1356	hydrodiag(:,:) = 0
1357	DO ig=1,nbpt
1358	DO ib=1,nbasmax
1359	ic = 0
1360	ign = ig
1361	ibn = ib
1362	! Locate outflow point
1363	DO WHILE (ibn .GT. 0 .AND. ibn .LE. nbasmax .AND. ic .LT. nbasmax*nbpt)
1364	ic = ic + 1
1365	og = ign
1366	ob = ibn
1367	ign = route_togrid(og, ob)
1368	ibn = route_tobasin(og, ob)
1369	ENDDO
1370	!
1371	! Now that we have an outflow check if it is one of the num_largest rivers.
1372	! In this case we keeps the location so we diagnose it.
1373	!
1374	IF ( ibn .EQ. nbasmax + 3) THEN
1375	DO icc = 1,num_largest
1376	IF ( outpt(icc,1) .EQ. og .AND. outpt(icc,2) .EQ. ob ) THEN
1377	!
1378	! We only keep this point for our map if it is large enough.
1379	!
1380	IF ( routing_area(ig,ib) .GT. 0.25resolution(ig,1)resolution(ig,2)*contfrac(ig) ) THEN
1381	nb_pts(icc) = nb_pts(icc) + 1
1382	pts(icc, nb_pts(icc)) = ig
1383	ptbas(icc, nb_pts(icc)) = ib
1384	ENDIF
1385	totarea(icc) = totarea(icc) + routing_area(ig,ib)
1386	! ID of the river is taken from the last point before the outflow.
1387	topids(icc) = global_basinid(og,ob)
1388	!
1389	! On this gridbox and basin we will diagnose the hydrograph
1390	!
1391	hydrodiag(ig, ib) = 1
1392	!
1393	ENDIF
1394	ENDDO
1395	ENDIF
1396
1397	ENDDO
1398	ENDDO
1399	!
1400	! Construct the map of the largest basins. We take the num_largest basins
1401	! if they have more than 4 points. After that it is of no use.
1402	!
1403	!
1404	basinmap(:) = zero
1405	DO icc = 1, num_largest
1406	ff = MAXLOC(totarea)
1407	IF ( nb_pts(ff(1)) .GT. 2 ) THEN
1408	DO ig = 1, nb_pts(ff(1))
1409	basinmap(pts(ff(1),ig)) = REAL(icc,r_std)
1410	ENDDO
1411	!
1412	IF ( topids(ff(1)) .GT. nb_small ) THEN
1413	WRITE(name_str, '("NN, Nb : ",I4)') topids(ff(1))
1414	ELSE
1415	name_str = basin_names(topids(ff(1)))
1416	ENDIF
1417	!
1418	WRITE(numout,&
1419	'("Basin ID ", I5," ", A15, " Area [km^2] : ", F13.4, " Nb points : ", I4)')&
1420	& topids(ff(1)), name_str(1:15), totarea(ff(1))/1.e6, nb_pts(ff(1))
1421	ENDIF
1422	totarea(ff(1)) = zero
1423	ENDDO
1424	!
1425	!
1426	nbrivers(:) = zero
1427	DO ig=1,nbpt
1428	nbrivers(ig) = COUNT(route_tobasin(ig,1:nbasmax) == nbasmax+3)
1429	ENDDO
1430	DO ig=1,nbpt
1431	IF ( nbrivers(ig) > 1 ) THEN
1432	WRITE(numout,*) 'Grid box ', ig, ' has ', NINT(nbrivers(ig)), ' outflow points.'
1433	WRITE(numout,*) 'The rivers which flow into the ocean at this point are :'
1434	DO icc=1,nbasmax
1435	IF ( route_tobasin(ig,icc) == nbasmax+3) THEN
1436	IF ( global_basinid(ig,icc) <= nb_small ) THEN
1437	WRITE(numout,*) 'ID = ',global_basinid(ig,icc), ' Name = ', basin_names(global_basinid(ig,icc))
1438	ELSE
1439	WRITE(numout,*) 'ID = ',global_basinid(ig,icc), ' Problem ===== ID is larger than possible'
1440	ENDIF
1441	ENDIF
1442	ENDDO
1443	ENDIF
1444	ENDDO
1445	!
1446	WRITE(numout,*) 'Maximum topographic index :', MAXVAL(topo_resid)
1447	ic = COUNT(topo_resid .GT. 0.)
1448	WRITE(numout,*) 'Mean topographic index :', SUM(topo_resid, MASK=topo_resid .GT. zero)/ic
1449	WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. zero)
1450	!
1451	DEALLOCATE(pts)
1452	DEALLOCATE(outpt)
1453	DEALLOCATE(nb_pts)
1454	DEALLOCATE(totarea)
1455	!
1456	END SUBROUTINE routing_diagnostic
1457	!
1458	!---------------------------------------------------------------------------------
1459	!
1460	SUBROUTINE routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac)
1461	!
1462	IMPLICIT NONE
1463	!
1464	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
1465	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
1466	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
1467	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
1468	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
1469
1470	! INTEGER(i_std) :: neighbours_tmp(nbpt,8)
1471	INTEGER(i_std) :: i,j
1472
1473	! DO i=1,nbp_loc
1474	! DO j=1,8
1475	! IF (neighbours(i,j)==-1) THEN
1476	! neighbours_tmp(i,j)=neighbours(i,j)
1477	! ELSE
1478	! neighbours_tmp(i,j)=neighbours(i,j)+nbp_para_begin(mpi_rank)-1
1479	! ENDIF
1480	! ENDDO
1481	! ENDDO
1482
1483	routing_area => routing_area_glo
1484	topo_resid => topo_resid_glo
1485	route_togrid => route_togrid_glo
1486	route_tobasin => route_tobasin_glo
1487	global_basinid => global_basinid_glo
1488
1489	IF (is_root_prc) CALL routing_basins(nbp_glo,lalo_g, neighbours_g, resolution_g, contfrac_g)
1490
1491	routing_area => routing_area_loc
1492	topo_resid => topo_resid_loc
1493	route_togrid => route_togrid_loc
1494	route_tobasin => route_tobasin_loc
1495	global_basinid => global_basinid_loc
1496
1497	CALL scatter(routing_area_glo,routing_area_loc)
1498	CALL scatter(topo_resid_glo,topo_resid_loc)
1499	CALL scatter(route_togrid_glo,route_togrid_loc)
1500	CALL scatter(route_tobasin_glo,route_tobasin_loc)
1501	CALL scatter(global_basinid_glo,global_basinid_loc)
1502
1503	END SUBROUTINE routing_basins_p
1504
1505	SUBROUTINE routing_basins(nbpt, lalo, neighbours, resolution, contfrac)
1506	!
1507	IMPLICIT NONE
1508	!
1509	!
1510	! This subroutine reads in the map of basins and flow direction to construct the
1511	! the catchments of each grid box.
1512	!
1513	! The work is done in a number of steps which are performed locally on the
1514	! GCM grid:
1515	! 1) First we find the grid-points of the high resolution routing grid which are
1516	! within the coarser grid of GCM.
1517	! 2) When we have these grid points we decompose them into basins in the routine
1518	! routing_findbasins. A number of simplifications are done if needed.
1519	! 3) In the routine routing_globalize we put the basin information of this grid
1520	! into global fields.
1521	! Then we work on the global grid to perform the following tasks :
1522	! 1) We linkup the basins of the various and check the global consistence.
1523	! 2) The area of each outflow point is computed.
1524	! 3) The final step is to reduce the number of basins in order to fit into the truncation.
1525	!
1526	!
1527	! 0.1 INPUT
1528	!
1529	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
1530	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
1531	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
1532	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
1533	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
1534	!
1535	!
1536	! 0.3 LOCAL
1537	!
1538	CHARACTER(LEN=80) :: filename
1539	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, lastjp, nbexp
1540	REAL(r_std) :: lev(1), date, dt, coslat
1541	INTEGER(i_std) :: itau(1), sgn
1542	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: trip, basins, topoindex, hierarchy
1543	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, lat_ful, lon_ful
1544	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: loup_rel, lolow_rel, laup_rel, lalow_rel
1545	REAL(r_std) :: lon_up, lon_low, lat_up, lat_low
1546	!
1547	INTEGER(i_std) :: nbi, nbj
1548	REAL(r_std) :: ax, ay, min_topoind, max_basins, invented_basins
1549	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
1550	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
1551	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: sub_pts
1552	REAL(r_std), DIMENSION(nbvmax,nbvmax) :: area_bx, hierarchy_bx, lon_bx, lat_bx, topoind_bx
1553	INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_bx, basin_bx
1554	!
1555	INTEGER(i_std) :: coast_pts(nbvmax)
1556	REAL(r_std) :: lonrel, louprel, lolowrel
1557	!
1558	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: basin_count
1559	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: basin_id
1560	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: basin_area, basin_hierarchy, basin_topoind
1561	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fetch_basin
1562	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: basin_flowdir
1563	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: outflow_grid
1564	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: outflow_basin
1565	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: inflow_number
1566	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: inflow_basin, inflow_grid
1567	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: nbcoastal
1568	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: coastal_basin
1569
1570	!
1571	INTEGER(i_std) :: nb_basin, nwbas
1572	INTEGER(i_std) :: basin_inbxid(nbvmax), basin_sz(nbvmax), basin_pts(nbvmax,nbvmax,2), basin_bxout(nbvmax)
1573	CHARACTER(LEN=7) :: fmt
1574	LOGICAL :: debug = .FALSE.
1575	INTEGER(i_std), DIMENSION(2) :: diagbox = (/ 1147, 1148 /)
1576
1577	! Test on diagbox and nbpt
1578	IF (debug) THEN
1579	IF (ANY(diagbox .GT. nbpt)) THEN
1580	WRITE(,) "Debug diganostics : nbpt, diagbox", nbpt, diagbox
1581	call ipslerr(3,'routing_basin', &
1582	& 'Problem with diagbox in debug mode.', &
1583	& 'diagbox values can''t be greater than land points number.', &
1584	& '(decrease diagbox wrong value)')
1585	ENDIF
1586	ENDIF
1587	!
1588	!
1589	! Needs to be a configurable variable
1590	!
1591	!
1592	!Config Key = ROUTING_FILE
1593	!Config Desc = Name of file which contains the routing information
1594	!Config Def = routing.nc
1595	!Config Help = The file provided here should alow the routing module to
1596	!Config read the high resolution grid of basins and the flow direction
1597	!Config from one mesh to the other.
1598	!
1599	filename = 'routing.nc'
1600	CALL getin('ROUTING_FILE',filename)
1601	!
1602	CALL flininfo(filename,iml, jml, lml, tml, fid)
1603	!
1604	!
1605	ALLOCATE (lat_rel(iml,jml))
1606	ALLOCATE (lon_rel(iml,jml))
1607	ALLOCATE (laup_rel(iml,jml))
1608	ALLOCATE (loup_rel(iml,jml))
1609	ALLOCATE (lalow_rel(iml,jml))
1610	ALLOCATE (lolow_rel(iml,jml))
1611	ALLOCATE (lat_ful(iml+2,jml+2))
1612	ALLOCATE (lon_ful(iml+2,jml+2))
1613	ALLOCATE (trip(iml,jml))
1614	ALLOCATE (basins(iml,jml))
1615	ALLOCATE (topoindex(iml,jml))
1616	ALLOCATE (hierarchy(iml,jml))
1617	!
1618	ALLOCATE (sub_area(nbpt,nbvmax))
1619	ALLOCATE (sub_index(nbpt,nbvmax,2))
1620	ALLOCATE (sub_pts(nbpt))
1621	!
1622
1623	!
1624	CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
1625	!
1626	! The trip field follows the following convention for the flow of the water :
1627	! trip = 1 : flow = N
1628	! trip = 2 : flow = NE
1629	! trip = 3 : flow = E
1630	! trip = 4 : flow = SE
1631	! trip = 5 : flow = S
1632	! trip = 6 : flow = SW
1633	! trip = 7 : flow = W
1634	! trip = 8 : flow = NW
1635	! trip = 97 : return flow into the ground
1636	! trip = 98 : coastal flow (diffuse flow into the oceans)
1637	! trip = 99 : river flow into the oceans
1638	!
1639	!
1640	CALL flinget(fid, 'trip', iml, jml, lml, tml, 1, 1, trip)
1641	!
1642	CALL flinget(fid, 'basins', iml, jml, lml, tml, 1, 1, basins)
1643	!
1644	CALL flinget(fid, 'topoind', iml, jml, lml, tml, 1, 1, topoindex)
1645	!
1646	CALL flinclo(fid)
1647	!
1648	nbexp = 0
1649	!
1650	min_topoind = MINVAL(topoindex, MASK=topoindex .LT. undef_sechiba-un)
1651	!
1652	DO ip=1,iml
1653	DO jp=1,jml
1654	IF ( trip(ip,jp) .LT. 1.e10 .AND. topoindex(ip,jp) .GT. 1.e10) THEN
1655	WRITE(numout,*) 'trip exists but not topoind :'
1656	WRITE(numout,*) 'ip, jp :', ip, jp
1657	WRITE(numout,*) 'trip, topoind : ', trip(ip,jp), topoindex(ip,jp)
1658	STOP
1659	ENDIF
1660	ENDDO
1661	ENDDO
1662	!
1663	! Duplicate the border assuming we have a global grid going from west to east
1664	!
1665	lon_ful(2:iml+1,2:jml+1) = lon_rel(1:iml,1:jml)
1666	lat_ful(2:iml+1,2:jml+1) = lat_rel(1:iml,1:jml)
1667	!
1668	IF ( lon_rel(iml,1) .LT. lon_ful(2,2)) THEN
1669	lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)
1670	lat_ful(1,2:jml+1) = lat_rel(iml,1:jml)
1671	ELSE
1672	lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)-360
1673	lat_ful(1,2:jml+1) = lat_rel(iml,1:jml)
1674	ENDIF
1675
1676	IF ( lon_rel(1,1) .GT. lon_ful(iml+1,2)) THEN
1677	lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)
1678	lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml)
1679	ELSE
1680	lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)+360
1681	lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml)
1682	ENDIF
1683	!
1684	sgn = INT(lat_rel(1,1)/ABS(lat_rel(1,1)))
1685	lat_ful(2:iml+1,1) = sgn*180 - lat_rel(1:iml,1)
1686	sgn = INT(lat_rel(1,jml)/ABS(lat_rel(1,jml)))
1687	lat_ful(2:iml+1,jml+2) = sgn*180 - lat_rel(1:iml,jml)
1688	lat_ful(1,1) = lat_ful(iml+1,1)
1689	lat_ful(iml+2,1) = lat_ful(2,1)
1690	lat_ful(1,jml+2) = lat_ful(iml+1,jml+2)
1691	lat_ful(iml+2,jml+2) = lat_ful(2,jml+2)
1692	!
1693	! Add the longitude lines to the top and bottom
1694	!
1695	lon_ful(:,1) = lon_ful(:,2)
1696	lon_ful(:,jml+2) = lon_ful(:,jml+1)
1697	!
1698	! Get the upper and lower limits of each grid box
1699	!
1700	DO ip=1,iml
1701	DO jp=1,jml
1702	!
1703	loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),&
1704	& 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
1705	lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),&
1706	& 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
1707	laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),&
1708	& 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
1709	lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),&
1710	& 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
1711	!
1712	ENDDO
1713	ENDDO
1714	!
1715	!
1716	!
1717	! Now we take each grid point and find out which values from the forcing we need to average
1718	!
1719	DO ib =1, nbpt
1720	!
1721	! We find the 4 limits of the grid-box. As we transform the resolution of the model
1722	! into longitudes and latitudes we do not have the problem of periodicity.
1723	! coslat is a help variable here !
1724	!
1725	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
1726	!
1727	lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat)
1728	lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
1729	!
1730	coslat = pi/180. * R_Earth
1731	!
1732	lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat)
1733	lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat)
1734	!
1735	!
1736	! Find the grid boxes from the data that go into the model's boxes
1737	! We still work as if we had a regular grid ! Well it needs to be localy regular so
1738	! so that the longitude at the latitude of the last found point is close to the one of the next point.
1739	!
1740	fopt = 0
1741	lastjp = 1
1742	DO ip=1,iml
1743	!
1744	! Either the center of the data grid point is in the interval of the model grid or
1745	! the East and West limits of the data grid point are on either sides of the border of
1746	! the data grid.
1747	!
1748	! To do that correctly we have to check if the grid box sits on the date-line.
1749	!
1750	IF ( lon_low < -180.0 ) THEN
1751	lonrel = MOD( lon_rel(ip,lastjp) - 360.0, 360.0)
1752	lolowrel = MOD( lolow_rel(ip,lastjp) - 360.0, 360.0)
1753	louprel = MOD( loup_rel(ip,lastjp) - 360.0, 360.0)
1754	!
1755	ELSE IF ( lon_up > 180.0 ) THEN
1756	lonrel = MOD( 360. - lon_rel(ip,lastjp), 360.0)
1757	lolowrel = MOD( 360. - lolow_rel(ip,lastjp), 360.0)
1758	louprel = MOD( 360. - loup_rel(ip,lastjp), 360.0)
1759	ELSE
1760	lonrel = lon_rel(ip,lastjp)
1761	lolowrel = lolow_rel(ip,lastjp)
1762	louprel = loup_rel(ip,lastjp)
1763	ENDIF
1764	!
1765	!
1766	!
1767	IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
1768	& lolowrel < lon_low .AND. louprel > lon_low .OR. &
1769	& lolowrel < lon_up .AND. louprel > lon_up ) THEN
1770	!
1771	DO jp = 1, jml
1772	!
1773	! Now that we have the longitude let us find the latitude
1774	!
1775	IF ( lat_rel(ip,jp) > lat_low .AND. lat_rel(ip,jp) < lat_up .OR. &
1776	& lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.&
1777	& lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN
1778	!
1779	lastjp = jp
1780	!
1781	! Is it a land point ?
1782	!
1783	IF (trip(ip,jp) .LT. 1.e10) THEN
1784	!
1785	fopt = fopt + 1
1786	IF ( fopt .GT. nbvmax) THEN
1787	WRITE(numout,*) 'Please increase nbvmax in subroutine routing_basins', ib
1788	STOP
1789	ELSE
1790	!
1791	! If we sit on the date line we need to do the same transformations as above.
1792	!
1793	IF ( lon_low < -180.0 ) THEN
1794	lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
1795	louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
1796	!
1797	ELSE IF ( lon_up > 180.0 ) THEN
1798	lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0)
1799	louprel = MOD( 360. - loup_rel(ip,jp), 360.0)
1800	ELSE
1801	lolowrel = lolow_rel(ip,jp)
1802	louprel = loup_rel(ip,jp)
1803	ENDIF
1804	!
1805	! Get the area of the fine grid in the model grid
1806	!
1807	coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), 0.001 )
1808	ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))pi/180. R_Earth * coslat
1809	ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))pi/180. R_Earth
1810	sub_area(ib, fopt) = ax*ay
1811	sub_index(ib, fopt, 1) = ip
1812	sub_index(ib, fopt, 2) = jp
1813	ENDIF
1814	ENDIF
1815	!
1816	sub_pts(ib) = fopt
1817	!
1818	ENDIF
1819	!
1820	ENDDO
1821	!
1822	ENDIF
1823	!
1824	ENDDO
1825	!
1826	!
1827	ENDDO
1828	!
1829	! Do some memory management.
1830	!
1831	DEALLOCATE (laup_rel)
1832	DEALLOCATE (loup_rel)
1833	DEALLOCATE (lalow_rel)
1834	DEALLOCATE (lolow_rel)
1835	DEALLOCATE (lat_ful)
1836	DEALLOCATE (lon_ful)
1837	!
1838	nwbas = MAXVAL(sub_pts)
1839	!
1840	ALLOCATE (basin_count(nbpt))
1841	ALLOCATE (basin_area(nbpt,nwbas), basin_hierarchy(nbpt,nwbas), basin_topoind(nbpt,nwbas))
1842	ALLOCATE (fetch_basin(nbpt,nwbas))
1843	ALLOCATE (basin_id(nbpt,nwbas), basin_flowdir(nbpt,nwbas))
1844	ALLOCATE (outflow_grid(nbpt,nwbas),outflow_basin(nbpt,nwbas))
1845	ALLOCATE (inflow_number(nbpt,nwbas))
1846	ALLOCATE (inflow_basin(nbpt,nwbas,nbvmax), inflow_grid(nbpt,nwbas,nbvmax))
1847	ALLOCATE (nbcoastal(nbpt), coastal_basin(nbpt,nwbas))
1848	!
1849	! Order all sub points in each grid_box and find the sub basins
1850	!
1851	! before we start we set the maps to empty
1852	!
1853	basin_id(:,:) = undef_int
1854	basin_count(:) = 0
1855	hierarchy(:,:) = undef_sechiba
1856	max_basins = MAXVAL(basins, MASK=basins .LT. 1.e10)
1857	invented_basins = max_basins
1858	nbcoastal(:) = 0
1859	!
1860	CALL routing_hierarchy(iml, jml, trip, topoindex, hierarchy)
1861	!
1862	!
1863	DO ib =1, nbpt
1864	!
1865	!
1866	! extract the information for this grid box
1867	!
1868	CALL routing_getgrid(nbpt, iml, jml, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
1869	& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, hierarchy, &
1870	& nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, hierarchy_bx, lon_bx, lat_bx)
1871	!
1872	CALL routing_findbasins(nbi, nbj, trip_bx, basin_bx, hierarchy_bx, topoind_bx,&
1873	& nb_basin, basin_inbxid, basin_sz, basin_bxout, basin_pts, coast_pts)
1874	!
1875	! Deal with the case where nb_basin=0 for this grid box. In this case all goes into coatal flow
1876	!
1877	IF ( debug .AND. (COUNT(diagbox .EQ. ib) .GT. 0) ) THEN
1878	WRITE(numout,*) '===================== IB = :', ib
1879	WRITE(numout,*) "sub_pts(ib) :", sub_pts(ib)
1880	WRITE(numout,*) 'LON LAT of GCM :', lalo(ib,2), lalo(ib,1)
1881	WRITE(numout,*) 'Neighbor options :', neighbours(ib,1:8)
1882	WRITE(numout,*) 'Resolution :', resolution(ib,1:2)
1883	WRITE(fmt,"('(',I3,'I6)')") nbi
1884	WRITE(numout,*) '-------------> trip ', trip_bx(1,1)
1885	DO jp=1,nbj
1886	WRITE(numout,fmt) trip_bx(1:nbi,jp)
1887	ENDDO
1888	WRITE(numout,*) '-------------> basin ',basin_bx(1,1)
1889	DO jp=1,nbj
1890	WRITE(numout,fmt) basin_bx(1:nbi,jp)
1891	ENDDO
1892	WRITE(numout,*) '-------------> hierarchy ',hierarchy_bx(1,1)
1893	DO jp=1,nbj
1894	WRITE(numout,fmt) INT(hierarchy_bx(1:nbi,jp)/mille)
1895	ENDDO
1896	WRITE(numout,*) '-------------> topoindex ',topoind_bx(1,1)
1897	DO jp=1,nbj
1898	WRITE(numout,fmt) INT(topoind_bx(1:nbi,jp)/mille)
1899	ENDDO
1900	!
1901	WRITE(numout,*) '------------> The basins we retain'
1902	DO jp=1,nb_basin
1903	WRITE(numout,*) 'index, size, bxout, coast :', basin_inbxid(jp), basin_sz(jp),&
1904	& basin_bxout(jp), coast_pts(jp)
1905	ENDDO
1906	!
1907	ENDIF
1908	!
1909	CALL routing_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarchy_bx, topoind_bx, min_topoind,&
1910	& nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,&
1911	& basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,&
1912	& nbcoastal, coastal_basin)
1913	!
1914	IF ( debug .AND. (COUNT(diagbox .EQ. ib) .GT. 0) ) THEN
1915	WRITE(numout,*) 'GLOBAL information after routing_globalize for grid ', ib
1916	DO jp=1,basin_count(ib)
1917	WRITE(numout,*) 'Basin ID : ', basin_id(ib, jp)
1918	WRITE(numout,*) 'Basin flowdir :', basin_flowdir(ib, jp)
1919	WRITE(numout,*) 'Basin hierarchy :', basin_hierarchy(ib, jp)
1920	WRITE(numout,*) 'Basin topoindex :', basin_topoind(ib, jp)
1921	WRITE(numout,*) 'Basin outflow grid :', outflow_grid(ib,jp)
1922	ENDDO
1923	ENDIF
1924	!
1925	ENDDO
1926	!
1927	CALL routing_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, &
1928	& basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, &
1929	& nbcoastal, coastal_basin, invented_basins)
1930	!
1931	WRITE(numout,*) 'The maximum number of basins in any grid :', MAXVAL(basin_count)
1932	!
1933	IF ( debug ) THEN
1934	DO ib=1,SIZE(diagbox)
1935	IF ( diagbox(ib) .GT. 0 ) THEN
1936	WRITE(numout,*) 'After routing_linkup information for grid ', diagbox(ib)
1937	DO jp=1,basin_count(diagbox(ib))
1938	WRITE(numout,*) 'Basin ID : ', basin_id(diagbox(ib), jp)
1939	WRITE(numout,*) 'Basin outflow_grid :', outflow_grid(diagbox(ib), jp)
1940	WRITE(numout,*) 'Basin outflow_basin:', outflow_basin(diagbox(ib), jp)
1941	WRITE(numout,*) 'Basin hierarchy :', basin_hierarchy(diagbox(ib), jp)
1942	ENDDO
1943	ENDIF
1944	ENDDO
1945	ENDIF
1946	!
1947	CALL routing_fetch(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_id, outflow_grid, &
1948	& outflow_basin, fetch_basin)
1949	!
1950	WRITE(numout,*) "Start reducing the number of basins per grid to meet the required truncation."
1951	!
1952	CALL routing_truncate(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_topoind,&
1953	& fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
1954	& inflow_grid, inflow_basin)
1955	!
1956	DEALLOCATE (lat_rel)
1957	DEALLOCATE (lon_rel)
1958	!
1959	DEALLOCATE (trip)
1960	DEALLOCATE (basins)
1961	DEALLOCATE (topoindex)
1962	DEALLOCATE (hierarchy)
1963	!
1964	DEALLOCATE (sub_area)
1965	DEALLOCATE (sub_index)
1966	DEALLOCATE (sub_pts)
1967	!
1968	DEALLOCATE (basin_count)
1969	DEALLOCATE (basin_area, basin_hierarchy, basin_topoind, fetch_basin)
1970	DEALLOCATE (basin_id, basin_flowdir)
1971	DEALLOCATE (outflow_grid,outflow_basin)
1972	DEALLOCATE (inflow_number)
1973	DEALLOCATE (inflow_basin, inflow_grid)
1974	DEALLOCATE (nbcoastal, coastal_basin)
1975	!
1976	RETURN
1977	!
1978	END SUBROUTINE routing_basins
1979	!
1980	!-----------------------------------------------------------------------
1981	!
1982	SUBROUTINE routing_getgrid(nbpt, iml, jml, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
1983	& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, hierarchy, &
1984	& nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, hierarchy_bx, lon_bx, lat_bx)
1985	!
1986	IMPLICIT NONE
1987	!
1988	! Extracts from the global high resolution fileds the data for the current grid box
1989	! we are dealing with.
1990	!
1991	! Convention for trip on the input :
1992	! The trip field follows the following convention for the flow of the water :
1993	! trip = 1 : flow = N
1994	! trip = 2 : flow = NE
1995	! trip = 3 : flow = E
1996	! trip = 4 : flow = SE
1997	! trip = 5 : flow = S
1998	! trip = 6 : flow = SW
1999	! trip = 7 : flow = W
2000	! trip = 8 : flow = NW
2001	! trip = 97 : return flow into the ground
2002	! trip = 98 : coastal flow (diffuse flow into the oceans) These values are created here
2003	! trip = 99 : river flow into the oceans
2004	!
2005	! On output, the grid boxes of the basin map which flow out of the GCM grid are identified
2006	! by numbers larger than 100 :
2007	! trip = 101 : flow = N out of the coarse grid
2008	! trip = 102 : flow = NE out of the coarse grid
2009	! trip = 103 : flow = E out of the coarse grid
2010	! trip = 104 : flow = SE out of the coarse grid
2011	! trip = 105 : flow = S out of the coarse grid
2012	! trip = 106 : flow = SW out of the coarse grid
2013	! trip = 107 : flow = W out of the coarse grid
2014	! trip = 108 : flow = NW out of the coarse grid
2015	! Inside the grid the convention remains the same as above (ie between 1 and 99).
2016	!
2017	! INPUT
2018	!
2019	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points in the global grid
2020	INTEGER(i_std), INTENT(in) :: iml, jml ! Resolution of the high resolution grid
2021	INTEGER(i_std), INTENT(in) :: ib ! point we are currently dealing with
2022	INTEGER(i_std), INTENT(in) :: sub_pts(nbpt) ! Number of high resiolution points on this grid
2023	INTEGER(i_std), INTENT(in) :: sub_index(nbpt, nbvmax,2) ! indeces of the points we need on the fine grid
2024	REAL(r_std), INTENT(inout) :: max_basins ! The current maximum of basins
2025	REAL(r_std), INTENT(in) :: min_topoind ! The current maximum of topographic index
2026	REAL(r_std), INTENT(in) :: sub_area(nbpt, nbvmax) ! area on the fine grid
2027	REAL(r_std), INTENT(in) :: lon_rel(iml, jml), lat_rel(iml, jml) ! coordinates of the fine grid
2028	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
2029	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y
2030	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box.
2031	REAL(r_std), INTENT(inout) :: trip(iml, jml), basins(iml, jml) ! data on the fine grid
2032	REAL(r_std), INTENT(inout) :: topoindex(iml, jml), hierarchy(iml, jml) ! data on the fine grid
2033	INTEGER(i_std), INTENT(out):: nbi, nbj ! Number of point ion x and y within the grid
2034	REAL(r_std), INTENT(out) :: area_bx(nbvmax,nbvmax), hierarchy_bx(nbvmax,nbvmax)
2035	REAL(r_std), INTENT(out) :: lon_bx(nbvmax,nbvmax), lat_bx(nbvmax,nbvmax), topoind_bx(nbvmax,nbvmax)
2036	INTEGER(i_std), INTENT(out):: trip_bx(nbvmax,nbvmax), basin_bx(nbvmax,nbvmax)
2037	!
2038	! LOCAL
2039	!
2040	INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc
2041	REAL(r_std) :: lonstr(nbvmaxnbvmax), latstr(nbvmaxnbvmax)
2042	!
2043	!
2044	! Set everything to undef to locate easily empty points
2045	!
2046	trip_bx(:,:) = undef_int
2047	basin_bx(:,:) = undef_int
2048	topoind_bx(:,:) = undef_sechiba
2049	area_bx(:,:) = undef_sechiba
2050	hierarchy_bx(:,:) = undef_sechiba
2051	!
2052	IF ( sub_pts(ib) > 0 ) THEN
2053	!
2054	DO ip=1,sub_pts(ib)
2055	lonstr(ip) = lon_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2))
2056	latstr(ip) = lat_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2))
2057	ENDDO
2058	!
2059	! Get the size of the area and order the coordinates to go from North to South and West to East
2060	!
2061	CALL routing_sortcoord(sub_pts(ib), lonstr, 'WE', nbi)
2062	CALL routing_sortcoord(sub_pts(ib), latstr, 'NS', nbj)
2063	!
2064	! Transfer the data in such a way that (1,1) is the North Western corner and
2065	! (nbi, nbj) the South Eastern.
2066	!
2067	DO ip=1,sub_pts(ib)
2068	ll = MINLOC(ABS(lonstr(1:nbi) - lon_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2))))
2069	iloc = ll(1)
2070	ll = MINLOC(ABS(latstr(1:nbj) - lat_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2))))
2071	jloc = ll(1)
2072	trip_bx(iloc, jloc) = NINT(trip(sub_index(ib, ip, 1), sub_index(ib, ip, 2)))
2073	basin_bx(iloc, jloc) = NINT(basins(sub_index(ib, ip, 1), sub_index(ib, ip, 2)))
2074	area_bx(iloc, jloc) = sub_area(ib, ip)
2075	topoind_bx(iloc, jloc) = topoindex(sub_index(ib, ip, 1), sub_index(ib, ip, 2))
2076	hierarchy_bx(iloc, jloc) = hierarchy(sub_index(ib, ip, 1), sub_index(ib, ip, 2))
2077	lon_bx(iloc, jloc) = lon_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2))
2078	lat_bx(iloc, jloc) = lat_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2))
2079	ENDDO
2080	ELSE
2081	!
2082	! This is the case where the model invented a continental point
2083	!
2084	nbi = 1
2085	nbj = 1
2086	iloc = 1
2087	jloc = 1
2088	trip_bx(iloc, jloc) = 98
2089	basin_bx(iloc, jloc) = NINT(max_basins + 1)
2090	max_basins = max_basins + 1
2091	area_bx(iloc, jloc) = resolution(ib,1)resolution(ib,2)contfrac(ib)
2092	topoind_bx(iloc, jloc) = min_topoind
2093	hierarchy_bx(iloc, jloc) = min_topoind
2094	lon_bx(iloc, jloc) = lalo(ib,2)
2095	lat_bx(iloc, jloc) = lalo(ib,1)
2096	!
2097	ENDIF
2098	!
2099	! Tag in trip all the outflow conditions. The table is thus :
2100	! trip = 100+n : Outflow into another grid-box
2101	! trip = 99 : River outflow into the ocean
2102	! trip = 98 : This will be coastal flow (not organized as a basin)
2103	! trip = 97 : return flow into the soil (local)
2104	!
2105	DO jp=1,nbj
2106	IF ( trip_bx(1,jp) .EQ. 8 .OR. trip_bx(1,jp) .EQ. 7 .OR. trip_bx(1,jp) .EQ. 6) THEN
2107	trip_bx(1,jp) = trip_bx(1,jp) + 100
2108	ENDIF
2109	IF ( trip_bx(nbi,jp) .EQ. 2 .OR. trip_bx(nbi,jp) .EQ. 3 .OR. trip_bx(nbi,jp) .EQ. 4) THEN
2110	trip_bx(nbi,jp) = trip_bx(nbi,jp) + 100
2111	ENDIF
2112	ENDDO
2113	DO ip=1,nbi
2114	IF ( trip_bx(ip,1) .EQ. 8 .OR. trip_bx(ip,1) .EQ. 1 .OR. trip_bx(ip,1) .EQ. 2) THEN
2115	trip_bx(ip,1) = trip_bx(ip,1) + 100
2116	ENDIF
2117	IF ( trip_bx(ip,nbj) .EQ. 6 .OR. trip_bx(ip,nbj) .EQ. 5 .OR. trip_bx(ip,nbj) .EQ. 4) THEN
2118	trip_bx(ip,nbj) = trip_bx(ip,nbj) + 100
2119	ENDIF
2120	ENDDO
2121	!
2122	!
2123	! We simplify the outflow. We only need the direction normal to the
2124	! box boundary and the 4 corners.
2125	!
2126	! Northern border
2127	IF ( trip_bx(1,1) .EQ. 102 ) trip_bx(1,1) = 101
2128	IF ( trip_bx(nbi,1) .EQ. 108 ) trip_bx(nbi,1) = 101
2129	DO ip=2,nbi-1
2130	IF ( trip_bx(ip,1) .EQ. 108 .OR. trip_bx(ip,1) .EQ. 102 ) trip_bx(ip,1) = 101
2131	ENDDO
2132	! Southern border
2133	IF ( trip_bx(1,nbj) .EQ. 104 ) trip_bx(1,nbj) = 105
2134	IF ( trip_bx(nbi,nbj) .EQ. 106 ) trip_bx(nbi,nbj) = 105
2135	DO ip=2,nbi-1
2136	IF ( trip_bx(ip,nbj) .EQ. 104 .OR. trip_bx(ip,nbj) .EQ. 106 ) trip_bx(ip,nbj) = 105
2137	ENDDO
2138	! Eastern border
2139	IF ( trip_bx(nbi,1) .EQ. 104) trip_bx(nbi,1) = 103
2140	IF ( trip_bx(nbi,nbj) .EQ. 102) trip_bx(nbi,nbj) = 103
2141	DO jp=2,nbj-1
2142	IF ( trip_bx(nbi,jp) .EQ. 104 .OR. trip_bx(nbi,jp) .EQ. 102 ) trip_bx(nbi,jp) = 103
2143	ENDDO
2144	! Western border
2145	IF ( trip_bx(1,1) .EQ. 106) trip_bx(1,1) = 107
2146	IF ( trip_bx(1,nbj) .EQ. 108) trip_bx(1,nbj) = 107
2147	DO jp=2,nbj-1
2148	IF ( trip_bx(1,jp) .EQ. 106 .OR. trip_bx(1,jp) .EQ. 108 ) trip_bx(1,jp) = 107
2149	ENDDO
2150	!
2151	!
2152	END SUBROUTINE routing_getgrid
2153	!
2154	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2155	!
2156	SUBROUTINE routing_sortcoord(nb_in, coords, direction, nb_out)
2157	!
2158	IMPLICIT NONE
2159	!
2160	! Input/Output
2161	!
2162	INTEGER(i_std), INTENT(in) :: nb_in
2163	REAL(r_std), INTENT(inout) :: coords(nb_in)
2164	CHARACTER(LEN=2) :: direction
2165	INTEGER(i_std), INTENT(out) :: nb_out
2166	!
2167	! Local
2168	!
2169	INTEGER(i_std) :: ipos
2170	REAL(r_std) :: coords_tmp(nb_in)
2171	INTEGER(i_std), DIMENSION(1) :: ll
2172	INTEGER(i_std) :: ind(nb_in)
2173	!
2174	ipos = 1
2175	nb_out = nb_in
2176	!
2177	! Compress the coordinates array
2178	!
2179	DO WHILE ( ipos < nb_in )
2180	IF ( coords(ipos+1) /= undef_sechiba) THEN
2181	IF ( COUNT(coords(ipos:nb_out) == coords(ipos)) > 1 ) THEN
2182	coords(ipos:nb_out-1) = coords(ipos+1:nb_out)
2183	coords(nb_out:nb_in) = undef_sechiba
2184	nb_out = nb_out - 1
2185	ELSE
2186	ipos = ipos + 1
2187	ENDIF
2188	ELSE
2189	EXIT
2190	ENDIF
2191	ENDDO
2192	!
2193	! Sort it now
2194	!
2195	! First we get ready and adjust for the periodicity in longitude
2196	!
2197	coords_tmp(:) = undef_sechiba
2198	IF ( INDEX(direction, 'WE') == 1 .OR. INDEX(direction, 'EW') == 1) THEN
2199	IF ( MAXVAL(ABS(coords(1:nb_out))) .GT. 160 ) THEN
2200	coords_tmp(1:nb_out) = MOD(coords(1:nb_out) + 360.0, 360.0)
2201	ELSE
2202	coords_tmp(1:nb_out) = coords(1:nb_out)
2203	ENDIF
2204	ELSE IF ( INDEX(direction, 'NS') == 1 .OR. INDEX(direction, 'SN') == 1) THEN
2205	coords_tmp(1:nb_out) = coords(1:nb_out)
2206	ELSE
2207	WRITE(numout,*) 'The chosen direction (', direction,') is not recognized'
2208	STOP 'routing_sortcoord'
2209	ENDIF
2210	!
2211	! Get it sorted out now
2212	!
2213	ipos = 1
2214	!
2215	IF ( INDEX(direction, 'WE') == 1 .OR. INDEX(direction, 'SN') == 1) THEN
2216	DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1)
2217	ll = MINLOC(coords_tmp(:), coords_tmp /= undef_sechiba)
2218	ind(ipos) = ll(1)
2219	coords_tmp(ll(1)) = undef_sechiba
2220	ipos = ipos + 1
2221	ENDDO
2222	ELSE IF ( INDEX(direction, 'EW') == 1 .OR. INDEX(direction, 'NS') == 1) THEN
2223	DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1)
2224	ll = MAXLOC(coords_tmp(:), coords_tmp /= undef_sechiba)
2225	ind(ipos) = ll(1)
2226	coords_tmp(ll(1)) = undef_sechiba
2227	ipos = ipos + 1
2228	ENDDO
2229	ELSE
2230	WRITE(numout,*) 'The chosen direction (', direction,') is not recognized (second)'
2231	STOP 'routing_sortcoord'
2232	ENDIF
2233	!
2234	coords(1:nb_out) = coords(ind(1:nb_out))
2235	IF (nb_out < nb_in) THEN
2236	coords(nb_out+1:nb_in) = zero
2237	ENDIF
2238	!
2239	END SUBROUTINE routing_sortcoord
2240	!
2241	!-------------------------------------------------------------------------------------------------
2242	!
2243	SUBROUTINE routing_findbasins(nbi, nbj, trip, basin, hierarchy, topoind, nb_basin, basin_inbxid, basin_sz,&
2244	& basin_bxout, basin_pts, coast_pts)
2245	!
2246	IMPLICIT NONE
2247	!
2248	! This subroutine find the basins and does some clean up. The aim is to return the list off all
2249	! points which are within the same basin of the grid_box.
2250	! We will also collect all points which directly flow into the ocean in one basin
2251	! Make sure that we do not have a basin with two outflows and other exceptions.
2252	!
2253	! At this stage no effort is made to come down to the truncation of the model.
2254	!
2255	! Convention for trip
2256	! -------------------
2257	! Inside of the box :
2258	! trip = 1 : flow = N
2259	! trip = 2 : flow = NE
2260	! trip = 3 : flow = E
2261	! trip = 4 : flow = SE
2262	! trip = 5 : flow = S
2263	! trip = 6 : flow = SW
2264	! trip = 7 : flow = W
2265	! trip = 8 : flow = NW
2266	! trip = 97 : return flow into the ground
2267	! trip = 98 : coastal flow (diffuse flow into the oceans) These values are created here
2268	! trip = 99 : river flow into the oceans
2269	!
2270	! Out flow from the gird :
2271	! trip = 101 : flow = N out of the coarse grid
2272	! trip = 102 : flow = NE out of the coarse grid
2273	! trip = 103 : flow = E out of the coarse grid
2274	! trip = 104 : flow = SE out of the coarse grid
2275	! trip = 105 : flow = S out of the coarse grid
2276	! trip = 106 : flow = SW out of the coarse grid
2277	! trip = 107 : flow = W out of the coarse grid
2278	! trip = 108 : flow = NW out of the coarse grid
2279	!
2280	! Inputs
2281	!
2282	INTEGER(i_std) :: nbi, nbj
2283	INTEGER(i_std) :: trip(:,:), basin(:,:)
2284	REAL(r_std) :: hierarchy(:,:), topoind(:,:)
2285	!
2286	! Outputs
2287	!
2288	INTEGER(i_std) :: nb_basin
2289	INTEGER(i_std) :: basin_inbxid(nbvmax), basin_sz(nbvmax), basin_bxout(nbvmax)
2290	INTEGER(i_std) :: basin_pts(nbvmax, nbvmax, 2)
2291	INTEGER(i_std) :: coast_pts(nbvmax)
2292	!
2293	! Local
2294	!
2295	INTEGER(i_std) :: ibas, ilf, nbb, nb_in
2296	INTEGER(i_std) :: bname(nbvmax), sz(nbvmax), pts(nbvmax,nbvmax,2), nbout(nbvmax)
2297	INTEGER(i_std) :: new_nb, new_bname(nbvmax), new_sz(nbvmax), new_pts(nbvmax,nbvmax,2)
2298	INTEGER(i_std) :: itrans, trans(nbvmax), outdir(nbvmax)
2299	INTEGER(i_std) :: tmpsz(nbvmax)
2300	INTEGER(i_std) :: ip, jp, jpp(1), ipb
2301	INTEGER(i_std) :: sortind(nbvmax)
2302	CHARACTER(LEN=7) :: fmt
2303	!
2304	nbb = 0
2305	bname(:) = undef_int
2306	sz(:) = 0
2307	nbout(:) = 0
2308	new_pts(:,:,:) = 0
2309	!
2310	! 1.0 Find all basins within this grid-box
2311	! Sort the variables per basin so that we can more easily
2312	! access data from the same basin (The variables are :
2313	! bname, sz, pts, nbout)
2314	!
2315	DO ip=1,nbi
2316	DO jp=1,nbj
2317	IF ( basin(ip,jp) .LT. undef_int) THEN
2318	IF ( COUNT(basin(ip,jp) .EQ. bname(:)) .EQ. 0 ) THEN
2319	nbb = nbb + 1
2320	IF ( nbb .GT. nbvmax ) STOP 'nbvmax too small'
2321	bname(nbb) = basin(ip,jp)
2322	sz(nbb) = 0
2323	ENDIF
2324	!
2325	DO ilf=1,nbb
2326	IF ( basin(ip,jp) .EQ. bname(ilf) ) THEN
2327	ibas = ilf
2328	ENDIF
2329	ENDDO
2330	!
2331	sz(ibas) = sz(ibas) + 1
2332	IF ( sz(ibas) .GT. nbvmax ) STOP 'nbvmax too small'
2333	pts(ibas, sz(ibas), 1) = ip
2334	pts(ibas, sz(ibas), 2) = jp
2335	! We deal only with outflow and leave flow back into the grid-box for later.
2336	IF ( trip(ip,jp) .GE. 97 ) THEN
2337	nbout(ibas) = nbout(ibas) + 1
2338	ENDIF
2339	!
2340	ENDIF
2341	!
2342	ENDDO
2343	ENDDO
2344	!
2345	! 2.0 All basins which have size 1 and flow to the ocean are put together.
2346	!
2347	itrans = 0
2348	coast_pts(:) = undef_int
2349	! Get all the points we can collect
2350	DO ip=1,nbb
2351	IF ( sz(ip) .EQ. 1 .AND. trip(pts(ip,1,1),pts(ip,1,2)) .EQ. 99) THEN
2352	itrans = itrans + 1
2353	trans(itrans) = ip
2354	trip(pts(ip,1,1),pts(ip,1,2)) = 98
2355	ENDIF
2356	ENDDO
2357	! put everything in the first basin
2358	IF ( itrans .GT. 1) THEN
2359	ipb = trans(1)
2360	coast_pts(sz(ipb)) = bname(ipb)
2361	bname(ipb) = -1
2362	DO ip=2,itrans
2363	sz(ipb) = sz(ipb) + 1
2364	coast_pts(sz(ipb)) = bname(trans(ip))
2365	sz(trans(ip)) = 0
2366	pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1)
2367	pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2)
2368	ENDDO
2369	ENDIF
2370	!
2371	! 3.0 Make sure that we have only one outflow point in each basin
2372	!
2373	! nbb is the number of basins on this grid box.
2374	new_nb = 0
2375	DO ip=1,nbb
2376	! We only do this for grid-points which have more than one outflow
2377	IF ( sz(ip) .GT. 1 .AND. nbout(ip) .GT. 1) THEN
2378	!
2379	! Pick up all points needed and store them in trans
2380	!
2381	itrans = 0
2382	DO jp=1,sz(ip)
2383	IF ( trip(pts(ip,jp,1),pts(ip,jp,2)) .GE. 97) THEN
2384	itrans = itrans + 1
2385	trans(itrans) = trip(pts(ip,jp,1),pts(ip,jp,2))
2386	ENDIF
2387	ENDDO
2388	!
2389	! First issue : We have more than one point of the basin which flows into
2390	! the ocean. In this case we put everything into coastal flow. It will go into
2391	! a separate basin in the routing_globalize routine.
2392	!
2393	IF ( (COUNT(trans(1:itrans) .EQ. 99) + COUNT(trans(1:itrans) .EQ. 98)) .GT. 1) THEN
2394	DO jp=1,sz(ip)
2395	IF ( trip(pts(ip,jp,1),pts(ip,jp,2)) .EQ. 99 ) THEN
2396	trip(pts(ip,jp,1),pts(ip,jp,2)) = 98
2397	trans(itrans) = trip(pts(ip,jp,1),pts(ip,jp,2))
2398	ENDIF
2399	ENDDO
2400	ENDIF
2401	!
2402	! Second issue : We have redundant outflows at the boundaries. That is two small grid
2403	! boxes flowing into the same GCM grid box.
2404	!
2405	IF ( COUNT(trans(1:itrans) .GT. 100) .GE. 1) THEN
2406	CALL routing_simplify(nbi, nbj, trip, basin, hierarchy, bname(ip))
2407	itrans = 0
2408	DO jp=1,sz(ip)
2409	IF ( trip(pts(ip,jp,1),pts(ip,jp,2)) .GE. 9) THEN
2410	itrans = itrans + 1
2411	trans(itrans) = trip(pts(ip,jp,1),pts(ip,jp,2))
2412	ENDIF
2413	ENDDO
2414	ENDIF
2415	!
2416	! Third issue : we have more than one outflow from the boxes. This could be
2417	! - flow into 2 or more neighboring GCM grids
2418	! - flow into a neighboring GCM grids and into the ocean or be a return flow (=97. =98, =99)
2419	! - flow into a neighboring GCM grids or ocean and back into the same GCM grid box
2420	! The only solution is to cut the basin up in as many parts.
2421	!
2422	IF ( COUNT(trans(1:itrans) .GE. 97) .GT. 1) THEN
2423	!
2424	nb_in = new_nb
2425	CALL routing_cutbasin(nbi, nbj, nbb, trip, basin, bname(ip), new_nb, new_bname, new_sz, new_pts)
2426	!
2427	! If we have split the basin then we need to cancel the old one
2428	!
2429	IF ( nb_in .NE. new_nb) THEN
2430	sz(ip) = 0
2431	ENDIF
2432	!
2433	ENDIF
2434	!
2435	ENDIF
2436	ENDDO
2437	!
2438	! Add the new basins to the end of the list
2439	!
2440	If ( nbb+new_nb .LE. nbvmax) THEN
2441	DO ip=1,new_nb
2442	bname(nbb+ip) = new_bname(ip)
2443	sz(nbb+ip) = new_sz(ip)
2444	pts(nbb+ip,:,:) = new_pts(ip,:,:)
2445	ENDDO
2446	nbb = nbb+new_nb
2447	ELSE
2448	WRITE(numout,*) 'Increase nbvmax. It is too small to contain all the basins (routing_findbasins)'
2449	STOP
2450	ENDIF
2451	!
2452	! Keep the output direction
2453	!
2454	DO ip=1,nbb
2455	IF ( sz(ip) .GT. 0 ) THEN
2456	trans(:) = 0
2457	DO jp=1,sz(ip)
2458	trans(jp) = trip(pts(ip,jp,1),pts(ip,jp,2))
2459	ENDDO
2460	outdir(ip) = MAXVAL(trans(1:sz(ip)))
2461	IF ( outdir(ip) .GE. 97 ) THEN
2462	outdir(ip) = outdir(ip) - 100
2463	ELSE
2464	WRITE(numout,*) 'Why are we here and can not find a trip larger than 96'
2465	WRITE(numout,*) 'Does this mean that the basin does not have any outflow ', ip, bname(ip)
2466	WRITE(fmt,"('(',I3,'I9)')") nbi
2467	WRITE(numout,*) '-----------------------> trip'
2468	DO jp=1,nbj
2469	WRITE(numout,fmt) trip(1:nbi,jp)
2470	ENDDO
2471	WRITE(numout,*) '-----------------------> basin'
2472	DO jp=1,nbj
2473	WRITE(numout,fmt) basin(1:nbi,jp)
2474	ENDDO
2475	STOP 'SUBROUTINE : routing_findbasins'
2476	ENDIF
2477	ENDIF
2478	ENDDO
2479	!
2480	!
2481	! Sort the output by size of the various basins.
2482	!
2483	nb_basin = COUNT(sz(1:nbb) .GT. 0)
2484	tmpsz(:) = -1
2485	tmpsz(1:nbb) = sz(1:nbb)
2486	DO ip=1,nbb
2487	jpp = MAXLOC(tmpsz(:))
2488	IF ( sz(jpp(1)) .GT. 0) THEN
2489	sortind(ip) = jpp(1)
2490	tmpsz(jpp(1)) = -1
2491	ENDIF
2492	ENDDO
2493	basin_inbxid(1:nb_basin) = bname(sortind(1:nb_basin))
2494	basin_sz(1:nb_basin) = sz(sortind(1:nb_basin))
2495	basin_pts(1:nb_basin,:,:) = pts(sortind(1:nb_basin),:,:)
2496	basin_bxout(1:nb_basin) = outdir(sortind(1:nb_basin))
2497	!
2498	! We can only check if we have at least as many outflows as basins
2499	!
2500	ip = COUNT(trip(1:nbi,1:nbj) .GE. 97 .AND. trip(1:nbi,1:nbj) .LT. undef_int)
2501	!! ip = ip + COUNT(trip(1:nbi,1:nbj) .EQ. 97)
2502	!! IF ( COUNT(trip(1:nbi,1:nbj) .EQ. 98) .GT. 0) ip = ip + 1
2503	IF ( ip .LT. nb_basin ) THEN
2504	WRITE(numout,*) 'We have less outflow points than basins :', ip
2505	WRITE(fmt,"('(',I3,'I9)')") nbi
2506	WRITE(numout,*) '-----------------------> trip'
2507	DO jp=1,nbj
2508	WRITE(numout,fmt) trip(1:nbi,jp)
2509	ENDDO
2510	WRITE(numout,*) '-----------------------> basin'
2511	DO jp=1,nbj
2512	WRITE(numout,fmt) basin(1:nbi,jp)
2513	ENDDO
2514	WRITE(numout,*) 'nb_basin :', nb_basin
2515	WRITE(numout,*) 'Basin sized :', basin_sz(1:nb_basin)
2516	STOP 'in routing_findbasins'
2517	ENDIF
2518	!
2519	END SUBROUTINE routing_findbasins
2520	!
2521	! ------------------------------------------------------------------------------------------
2522	!
2523	SUBROUTINE routing_simplify(nbi, nbj, trip, basin, hierarchy, basin_inbxid)
2524	!
2525	IMPLICIT NONE
2526	!
2527	! This subroutine symplifies the routing out of each basin by taking
2528	! out redundancies at the borders of the GCM box. The aim is to have
2529	! only one outflow point per basin. But here we will not change the
2530	! the direction of the outflow.
2531	!
2532	! Inputs
2533	INTEGER(i_std) :: nbi, nbj
2534	INTEGER(i_std) :: trip(:,:), basin(:,:)
2535	REAL(r_std) :: hierarchy(:,:)
2536	INTEGER(i_std) :: basin_inbxid
2537	!
2538	! Local
2539	!
2540	INTEGER(i_std) :: ip, jp, nbout, basin_sz, iborder
2541	INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_tmp
2542	INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: trip_flow
2543	INTEGER(i_std), DIMENSION(nbvmax,2) :: outflow
2544	INTEGER(i_std), DIMENSION(nbvmax) :: outsz
2545	CHARACTER(LEN=7) :: fmt
2546	!
2547	INTEGER(i_std), DIMENSION(8,2) :: inc
2548	INTEGER(i_std) :: itodo, ill(1), icc, ismall, ibas, iip, jjp, ib, id
2549	INTEGER(i_std), DIMENSION(nbvmax) :: todopt, todosz
2550	REAL(r_std), DIMENSION(nbvmax) :: todohi
2551	LOGICAL :: not_found, debug = .FALSE.
2552	!
2553	!
2554	! The routing code (i=1, j=2)
2555	!
2556	inc(1,1) = 0
2557	inc(1,2) = -1
2558	inc(2,1) = 1
2559	inc(2,2) = -1
2560	inc(3,1) = 1
2561	inc(3,2) = 0
2562	inc(4,1) = 1
2563	inc(4,2) = 1
2564	inc(5,1) = 0
2565	inc(5,2) = 1
2566	inc(6,1) = -1
2567	inc(6,2) = 1
2568	inc(7,1) = -1
2569	inc(7,2) = 0
2570	inc(8,1) = -1
2571	inc(8,2) = -1
2572	!
2573	!
2574	! Symplify the outflow conditions first. We are only interested in the
2575	! outflows which go to different GCM grid-boxes.
2576	!
2577	IF ( debug ) THEN
2578	WRITE(numout,*) '+++++++++++++++++++ BEFORE ANYTHING ++++++++++++++++++++'
2579	WRITE(fmt,"('(',I3,'I6)')") nbi
2580	DO jp=1,nbj
2581	WRITE(numout,fmt) trip_tmp(1:nbi,jp)
2582	ENDDO
2583	ENDIF
2584	!
2585	! transfer the trips into an array which only contains the basin we are interested in
2586	!
2587	trip_tmp(:,:) = -1
2588	basin_sz = 0
2589	DO ip=1,nbi
2590	DO jp=1,nbj
2591	IF ( basin(ip,jp) .EQ. basin_inbxid) THEN
2592	trip_tmp(ip,jp) = trip(ip,jp)
2593	basin_sz = basin_sz + 1
2594	ENDIF
2595	ENDDO
2596	ENDDO
2597	!
2598	! Determine for each point where it flows to
2599	!
2600	CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz)
2601	!
2602	!
2603	!
2604	!
2605	! Over the width of a GCM grid box we can have many outflows but we are interested
2606	! in only one for each basin. Thus we wish to collect them all to form only one outflow
2607	! to the neighboring grid-box.
2608	!
2609	DO iborder = 101,107,2
2610	!
2611	! If we have more than one of these outflows then we need to merge the sub-basins
2612	!
2613	icc = COUNT(trip_tmp .EQ. iborder)-1
2614	DO WHILE ( icc .GT. 0)
2615	! Pick out all the points we will have to do
2616	itodo = 0
2617	DO ip=1,nbout
2618	IF (trip_tmp(outflow(ip,1),outflow(ip,2)) .EQ. iborder) THEN
2619	itodo = itodo + 1
2620	todopt(itodo) = ip
2621	todosz(itodo) = outsz(ip)
2622	! We take the hierarchy of the outflow point as we will try to
2623	! minimize if for the outflow of the entire basin.
2624	todohi(itodo) = hierarchy(outflow(ip,1),outflow(ip,2))
2625	ENDIF
2626	ENDDO
2627	!
2628	! We change the direction of the smalest basin.
2629	!
2630	ill=MAXLOC(todohi(1:itodo))
2631	ismall = todopt(ill(1))
2632	!
2633	DO ip=1,nbi
2634	DO jp=1,nbj
2635	IF ( trip_flow(ip,jp,1) .EQ. outflow(ismall,1) .AND.&
2636	& trip_flow(ip,jp,2) .EQ. outflow(ismall,2) ) THEN
2637	! Now that we have found a point of the smallest sub-basin we
2638	! look around for another sub-basin
2639	ib = 1
2640	not_found = .TRUE.
2641	DO WHILE ( not_found .AND. ib .LE. itodo )
2642	IF ( ib .NE. ill(1) ) THEN
2643	ibas = todopt(ib)
2644	DO id=1,8
2645	iip = ip + inc(id,1)
2646	jjp = jp + inc(id,2)
2647	! Can we look at this points or is there any need to ?
2648	IF ( iip .GE. 1 .AND. iip .LE. nbi .AND. &
2649	& jjp .GE. 1 .AND. jjp .LE. nbj .AND. not_found) THEN
2650	! Is this point the one we look for ?
2651	IF ( trip_flow(iip,jjp,1) .EQ. outflow(ibas,1) .AND. &
2652	& trip_flow(iip,jjp,2) .EQ. outflow(ibas,2)) THEN
2653	trip_flow(ip,jp,1) = outflow(ibas,1)
2654	trip_flow(ip,jp,2) = outflow(ibas,2)
2655	trip_tmp(ip,jp) = id
2656	! This last line ensures that we do not come back to this point
2657	! and that in the end the outer while will stop
2658	not_found = .FALSE.
2659	ENDIF
2660	ENDIF
2661	ENDDO
2662	ENDIF
2663	ib = ib + 1
2664	ENDDO
2665	ENDIF
2666	ENDDO
2667	ENDDO
2668	!
2669	icc = icc - 1
2670	ENDDO
2671	!
2672	!
2673	ENDDO
2674	!
2675	IF ( debug ) THEN
2676	WRITE(numout,*) '+++++++++++++++++++ AFTER +++++++++++++++++++++++++++++'
2677	WRITE(fmt,"('(',I3,'I6)')") nbi
2678	DO jp=1,nbj
2679	WRITE(numout,fmt) trip_tmp(1:nbi,jp)
2680	ENDDO
2681	ENDIF
2682	!
2683	! Put trip_tmp back into trip
2684	!
2685	DO ip=1,nbi
2686	DO jp=1,nbj
2687	IF ( trip_tmp(ip,jp) .GT. 0) THEN
2688	trip(ip,jp) = trip_tmp(ip,jp)
2689	ENDIF
2690	ENDDO
2691	ENDDO
2692	!
2693	END SUBROUTINE routing_simplify
2694	!
2695	!-------------------------------------
2696	!
2697	SUBROUTINE routing_cutbasin (nbi, nbj, nbbasins, trip, basin, basin_inbxid, nb, bname, sz, pts)
2698	!
2699	IMPLICIT NONE
2700	!
2701	! This subroutine cuts the original basin which has more than one outflow into as
2702	! many subbasins as outflow directions.
2703	!
2704	! Inputs
2705	INTEGER(i_std) :: nbi, nbj
2706	INTEGER(i_std) :: nbbasins
2707	INTEGER(i_std) :: trip(:,:), basin(:,:)
2708	INTEGER(i_std) :: basin_inbxid
2709	!
2710	! Outputs
2711	!
2712	INTEGER(i_std) :: nb, bname(nbvmax), sz(nbvmax), pts(nbvmax,nbvmax,2)
2713	!
2714	! Local
2715	!
2716	INTEGER(i_std) :: ip, jp, iip, jjp, ib, ibb, id, nbout
2717	INTEGER(i_std) :: basin_sz, nb_in
2718	INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_tmp
2719	INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: trip_flow
2720	INTEGER(i_std), DIMENSION(nbvmax,2) :: outflow
2721	INTEGER(i_std), DIMENSION(nbvmax) :: outsz
2722	CHARACTER(LEN=7) :: fmt
2723	LOGICAL :: not_found, debug=.FALSE.
2724	!
2725	INTEGER(i_std), DIMENSION(8,2) :: inc
2726	!
2727	!
2728	! The routing code (i=1, j=2)
2729	!
2730	inc(1,1) = 0
2731	inc(1,2) = -1
2732	inc(2,1) = 1
2733	inc(2,2) = -1
2734	inc(3,1) = 1
2735	inc(3,2) = 0
2736	inc(4,1) = 1
2737	inc(4,2) = 1
2738	inc(5,1) = 0
2739	inc(5,2) = 1
2740	inc(6,1) = -1
2741	inc(6,2) = 1
2742	inc(7,1) = -1
2743	inc(7,2) = 0
2744	inc(8,1) = -1
2745	inc(8,2) = -1
2746	!
2747	! Set up a temporary trip field which only contains the values
2748	! for the basin on which we currently work.
2749	!
2750	trip_tmp(:,:) = -1
2751	basin_sz = 0
2752	DO ip=1,nbi
2753	DO jp=1,nbj
2754	IF ( basin(ip,jp) .EQ. basin_inbxid) THEN
2755	trip_tmp(ip,jp) = trip(ip,jp)
2756	basin_sz = basin_sz + 1
2757	ENDIF
2758	ENDDO
2759	ENDDO
2760	!
2761	CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz)
2762	!
2763	!! IF ( debug ) THEN
2764	!! DO ib = nb_in+1,nb
2765	!! DO ip=1,sz(ib)
2766	!! trip_tmp(pts(ib, ip, 1),pts(ib, ip, 2)) = ib*(-1)-900
2767	!! ENDDO
2768	!! ENDDO
2769	!! WRITE(fmt,"('(',I3,'I6)')") nbi
2770	!! WRITE(numout,*) 'BEFORE ------------> New basins '
2771	!! WRITE(numout,*) nb, ' sz :', sz(1:nb)
2772	!! DO jp=1,nbj
2773	!! WRITE(numout,fmt) trip_tmp(1:nbi,jp)
2774	!! ENDDO
2775	!! ENDIF
2776	!
2777	! Take out the small sub-basins. That is those which have only one grid box
2778	! This is only done if we need to save space in the number of basins. Else we
2779	! can take it easy and keep diverging sub-basins for the moment.
2780	!
2781	IF ( nbbasins .GE. nbasmax ) THEN
2782	DO ib=1,nbout
2783	! If the sub-basin is of size one and its larger neighbor is flowing into another
2784	! direction then we put them together.
2785	IF ( outsz(ib) .EQ. 1 .AND. trip(outflow(ib,1), outflow(ib,2)) .GT. 99 ) THEN
2786	!
2787	not_found = .TRUE.
2788	DO id=1,8
2789	ip = outflow(ib,1)
2790	jp = outflow(ib,2)
2791	iip = ip + inc(id,1)
2792	jjp = jp + inc(id,2)
2793	! Can we look at this points ?
2794	IF ( iip .GE. 1 .AND. iip .LE. nbi .AND. &
2795	& jjp .GE. 1 .AND. jjp .LE. nbj .AND. not_found) THEN
2796	! Did we find a direct neighbor which is an outflow point ?
2797	IF ( trip_tmp(iip,jjp) .GT. 100 ) THEN
2798	! IF so direct the flow towards it and update the tables.
2799	not_found = .FALSE.
2800	trip(ip,jp) = id
2801	trip_tmp(ip,jp) = id
2802	outsz(ib) = 0
2803	! update the table of this basin
2804	DO ibb=1,nbout
2805	IF ( iip .EQ. outflow(ibb,1) .AND. jjp .EQ. outflow(ibb,2) ) THEN
2806	outsz(ibb) = outsz(ibb)+1
2807	trip_flow(ip,jp,1) = outflow(ibb,1)
2808	trip_flow(ip,jp,2) = outflow(ibb,2)
2809	ENDIF
2810	ENDDO
2811	ENDIF
2812	ENDIF
2813	ENDDO
2814	ENDIF
2815	ENDDO
2816	ENDIF
2817	!
2818	!
2819	! Cut the basin if we have more than 1 left.
2820	!
2821	!
2822	IF ( COUNT(outsz(1:nbout) .GT. 0) .GT. 1 ) THEN
2823	!
2824	nb_in = nb
2825	!
2826	DO ib = 1,nbout
2827	IF ( outsz(ib) .GT. 0) THEN
2828	nb = nb+1
2829	IF ( nb .GT. nbvmax) THEN
2830	WRITE(numout,*) 'nbvmax too small, increase it (routing_cutbasin)'
2831	ENDIF
2832	bname(nb) = basin_inbxid
2833	sz(nb) = 0
2834	DO ip=1,nbi
2835	DO jp=1,nbj
2836	IF ( (trip_flow(ip,jp,1) + trip_flow(ip,jp,1)) .GT. 0 .AND. &
2837	& trip_flow(ip,jp,1) .EQ. outflow(ib,1) .AND. &
2838	& trip_flow(ip,jp,2) .EQ. outflow(ib,2) ) THEN
2839	sz(nb) = sz(nb) + 1
2840	pts(nb, sz(nb), 1) = ip
2841	pts(nb, sz(nb), 2) = jp
2842	ENDIF
2843	ENDDO
2844	ENDDO
2845	ENDIF
2846	ENDDO
2847	! A short verification
2848	IF ( SUM(sz(nb_in+1:nb)) .NE. basin_sz) THEN
2849	WRITE(numout,*) 'Lost some points while spliting the basin'
2850	WRITE(numout,*) 'nbout :', nbout
2851	DO ib = nb_in+1,nb
2852	WRITE(numout,*) 'ib, SZ :', ib, sz(ib)
2853	ENDDO
2854	WRITE(fmt,"('(',I3,'I6)')") nbi
2855	WRITE(numout,*) '-------------> trip '
2856	DO jp=1,nbj
2857	WRITE(numout,fmt) trip_tmp(1:nbi,jp)
2858	ENDDO
2859	STOP
2860	ENDIF
2861	!
2862	IF ( debug ) THEN
2863	DO ib = nb_in+1,nb
2864	DO ip=1,sz(ib)
2865	trip_tmp(pts(ib, ip, 1),pts(ib, ip, 2)) = ib*(-1)-900
2866	ENDDO
2867	ENDDO
2868	WRITE(fmt,"('(',I3,'I6)')") nbi
2869	WRITE(numout,*) 'AFTER-------------> New basins '
2870	WRITE(numout,*) nb, ' sz :', sz(1:nb)
2871	DO jp=1,nbj
2872	WRITE(numout,fmt) trip_tmp(1:nbi,jp)
2873	ENDDO
2874	IF ( MAXVAl(trip_tmp(1:nbi,1:nbj)) .GT. 0) THEN
2875	STOP
2876	ENDIF
2877	ENDIF
2878	ENDIF
2879	!
2880	END SUBROUTINE routing_cutbasin
2881	!
2882	!------------------------------------------------------------------------
2883	!
2884	SUBROUTINE routing_hierarchy(iml, jml, trip, topoindex, hierarchy)
2885	!
2886	IMPLICIT NONE
2887	!
2888	! This subroutine will find for each point the distance to the outflow point
2889	! along the flowlines of the basin.
2890	!
2891	INTEGER(i_std) :: iml, jml
2892	REAL(r_std), DIMENSION(iml,jml) :: trip, hierarchy, topoindex
2893	!
2894	INTEGER(i_std), DIMENSION(8,2) :: inc
2895	INTEGER(i_std) :: ip, jp, ib, ntripi, ntripj, cnt, trp
2896	REAL(r_std) :: topohier, topohier_old
2897	CHARACTER(LEN=7) :: fmt
2898	!
2899	! The routing code (i=1, j=2)
2900	!
2901	inc(1,1) = 0
2902	inc(1,2) = -1
2903	inc(2,1) = 1
2904	inc(2,2) = -1
2905	inc(3,1) = 1
2906	inc(3,2) = 0
2907	inc(4,1) = 1
2908	inc(4,2) = 1
2909	inc(5,1) = 0
2910	inc(5,2) = 1
2911	inc(6,1) = -1
2912	inc(6,2) = 1
2913	inc(7,1) = -1
2914	inc(7,2) = 0
2915	inc(8,1) = -1
2916	inc(8,2) = -1
2917	!
2918	DO ip=1,iml
2919	DO jp=1,jml
2920	IF ( trip(ip,jp) .LT. undef_sechiba ) THEN
2921	ntripi = ip
2922	ntripj = jp
2923	trp = trip(ip,jp)
2924	cnt = 1
2925	! Warn for extreme numbers
2926	IF ( topoindex(ip,jp) .GT. 1.e10 ) THEN
2927	WRITE(numout,*) 'We have a very large topographic index for point ', ip, jp
2928	WRITE(numout,*) 'This can not be right :', topoindex(ip,jp)
2929	STOP
2930	ELSE
2931	topohier = topoindex(ip,jp)
2932	ENDIF
2933	!
2934	DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml)
2935	cnt = cnt + 1
2936	ntripi = ntripi + inc(trp,1)
2937	IF ( ntripi .LT. 1) ntripi = iml
2938	IF ( ntripi .GT. iml) ntripi = 1
2939	ntripj = ntripj + inc(trp,2)
2940	topohier_old = topohier
2941	topohier = topohier + topoindex(ntripi, ntripj)
2942	IF ( topohier_old .GT. topohier) THEN
2943	WRITE(numout,*) 'Big Problem, how comes we climb up a hill ?'
2944	WRITE(numout,*) 'The old value of topographicaly weighted hierarchy was : ', topohier_old
2945	WRITE(numout,*) 'The new one is :', topohier
2946	STOP 'routing_hierarchy'
2947	ENDIF
2948	trp = trip(ntripi, ntripj)
2949	ENDDO
2950	!
2951	IF ( cnt .EQ. iml*jml) THEN
2952	WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp
2953	WRITE(numout,*) '-------------> trip '
2954	WRITE(fmt,"('(',I3,'I6)')") iml
2955	DO ib=1,jml
2956	WRITE(numout,fmt) trip(1:iml,ib)
2957	ENDDO
2958	STOP
2959	ENDIF
2960	!
2961	hierarchy(ip,jp) = topohier
2962	!
2963	ENDIF
2964	ENDDO
2965	ENDDO
2966	!
2967	!
2968	END SUBROUTINE routing_hierarchy
2969	!
2970	!------------------------------------------------------------------------
2971	!
2972	SUBROUTINE routing_findrout(nbi, nbj, trip, basin_sz, basinid, nbout, outflow, trip_flow, outsz)
2973	!
2974	IMPLICIT NONE
2975	!
2976	! This subroutine simply computes the route to each outflow point and returns
2977	! the outflow point for each point in the basin.
2978	!
2979	! INPUT
2980	!
2981	INTEGER(i_std) :: nbi, nbj
2982	INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip
2983	INTEGER(i_std) :: basin_sz, basinid, nbout
2984	INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: trip_flow
2985	INTEGER(i_std), DIMENSION(nbvmax,2) :: outflow
2986	INTEGER(i_std), DIMENSION(nbvmax) :: outsz
2987	!
2988	! LOCAL
2989	!
2990	INTEGER(i_std), DIMENSION(8,2) :: inc
2991	INTEGER(i_std) :: ip, jp, ib, cnt, trp, totsz
2992	CHARACTER(LEN=7) :: fmt
2993	!
2994	!
2995	! The routing code (i=1, j=2)
2996	!
2997	inc(1,1) = 0
2998	inc(1,2) = -1
2999	inc(2,1) = 1
3000	inc(2,2) = -1
3001	inc(3,1) = 1
3002	inc(3,2) = 0
3003	inc(4,1) = 1
3004	inc(4,2) = 1
3005	inc(5,1) = 0
3006	inc(5,2) = 1
3007	inc(6,1) = -1
3008	inc(6,2) = 1
3009	inc(7,1) = -1
3010	inc(7,2) = 0
3011	inc(8,1) = -1
3012	inc(8,2) = -1
3013	!
3014	!
3015	! Get the outflows and determine for each point to which outflow point it belong
3016	!
3017	nbout = 0
3018	trip_flow(:,:,:) = 0
3019	DO ip=1,nbi
3020	DO jp=1,nbj
3021	IF ( trip(ip,jp) .GT. 9) THEN
3022	nbout = nbout + 1
3023	outflow(nbout,1) = ip
3024	outflow(nbout,2) = jp
3025	ENDIF
3026	IF ( trip(ip,jp) .GT. 0) THEN
3027	trip_flow(ip,jp,1) = ip
3028	trip_flow(ip,jp,2) = jp
3029	ENDIF
3030	ENDDO
3031	ENDDO
3032	!
3033	! Follow the flow of the water
3034	!
3035	DO ip=1,nbi
3036	DO jp=1,nbj
3037	IF ( (trip_flow(ip,jp,1) + trip_flow(ip,jp,2)) .GT. 0) THEN
3038	trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2))
3039	cnt = 0
3040	DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj)
3041	cnt = cnt + 1
3042	trip_flow(ip,jp,1) = trip_flow(ip,jp,1) + inc(trp,1)
3043	trip_flow(ip,jp,2) = trip_flow(ip,jp,2) + inc(trp,2)
3044	trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2))
3045	ENDDO
3046	IF ( cnt .EQ. nbi*nbj) THEN
3047	WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp
3048	WRITE(numout,*) '-------------> trip '
3049	WRITE(fmt,"('(',I3,'I6)')") nbi
3050	DO ib=1,nbj
3051	WRITE(numout,fmt) trip(1:nbi,ib)
3052	ENDDO
3053	STOP
3054	ENDIF
3055	ENDIF
3056	ENDDO
3057	ENDDO
3058	!
3059	! What is the size of the region behind each outflow point ?
3060	!
3061	totsz = 0
3062	DO ip=1,nbout
3063	outsz(ip) = COUNT(trip_flow(:,:,1) .EQ. outflow(ip,1) .AND. trip_flow(:,:,2) .EQ. outflow(ip,2))
3064	totsz = totsz + outsz(ip)
3065	ENDDO
3066	IF ( basin_sz .NE. totsz) THEN
3067	WRITE(numout,*) 'Water got lost while I tried to follow it '
3068	WRITE(numout,*) basin_sz, totsz
3069	WRITE(numout,*) 'Basin id :', basinid
3070	DO ip=1,nbout
3071	WRITE(numout,*) 'ip :', ip, ' outsz :', outsz(ip), ' outflow :', outflow(ip,1), outflow(ip,2)
3072	ENDDO
3073	WRITE(numout,*) '-------------> trip '
3074	WRITE(fmt,"('(',I3,'I6)')") nbi
3075	DO jp=1,nbj
3076	WRITE(numout,fmt) trip(1:nbi,jp)
3077	ENDDO
3078	STOP
3079	ENDIF
3080	!
3081	END SUBROUTINE routing_findrout
3082	!
3083	!----------------------------------------------------------------------------------------------
3084	!
3085	SUBROUTINE routing_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarchy_bx, topoind_bx, min_topoind,&
3086	& nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,&
3087	& basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,&
3088	& nbcoastal, coastal_basin)
3089	!
3090	IMPLICIT NONE
3091	!
3092	! This subroutine will put the basins found for grid-box in in the global map. Connection can
3093	! only be made later when all information is together.
3094	!
3095	!
3096	! One of the outputs is basin_flowdir. Its convention is 1-8 for the directions from North to North
3097	! West going through South. The negative values will be -3 for return flow, -2 for coastal flow
3098	! and -1 river flow.
3099	!
3100	! LOCAL
3101	!
3102	INTEGER(i_std), INTENT (in) :: nbpt, ib ! Current grid-box
3103	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point
3104	! (1=N, 2=E, 3=S, 4=W)
3105	REAL(r_std), DIMENSION(nbvmax,nbvmax) :: area_bx ! Area of each small box in the grid-box
3106	INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_bx ! The trip field for each of the smaler boxes
3107	REAL(r_std), DIMENSION(nbvmax,nbvmax) :: hierarchy_bx ! level in the basin of the point
3108	REAL(r_std), DIMENSION(nbvmax,nbvmax) :: topoind_bx ! Topographic index
3109	REAL(r_std) :: min_topoind ! The default topographic index
3110	INTEGER(i_std) :: nb_basin ! number of sub-basins
3111	INTEGER(i_std), DIMENSION(nbvmax) :: basin_inbxid, basin_sz ! ID of basin, number of points in the basin
3112	INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2):: basin_pts ! Points in each basin
3113	INTEGER(i_std), DIMENSION(nbvmax) :: basin_bxout ! outflow direction
3114	INTEGER(i_std) :: coast_pts(nbvmax) ! The coastal flow points
3115	! global maps
3116	INTEGER(i_std) :: nwbas
3117	INTEGER(i_std), DIMENSION(nbpt) :: basin_count
3118	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id, basin_flowdir
3119	REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area, basin_hierarchy, basin_topoind
3120	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid
3121	INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal
3122	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin
3123	!
3124	! LOCAL
3125	!
3126	INTEGER(i_std) :: ij, iz
3127	CHARACTER(LEN=4) :: hierar_method = 'OUTP'
3128	!
3129	!
3130	DO ij=1, nb_basin
3131	!
3132	! Count the basins and keep their ID
3133	!
3134	basin_count(ib) = basin_count(ib)+1
3135	if (basin_count(ib) > nwbas) then
3136	WRITE(numout,*) 'ib=',ib
3137	call ipslerr(3,'routing_globalize', &
3138	& 'Problem with basin_count : ', &
3139	& 'It is greater than number of allocated basin nwbas.', &
3140	& '(stop to count basins)')
3141	endif
3142	basin_id(ib,basin_count(ib)) = basin_inbxid(ij)
3143	!
3144	! Transfer the list of basins which flow into the ocean as coastal flow.
3145	!
3146	IF ( basin_id(ib,basin_count(ib)) .LT. 0) THEN
3147	nbcoastal(ib) = basin_sz(ij)
3148	coastal_basin(ib,1:nbcoastal(ib)) = coast_pts(1:nbcoastal(ib))
3149	ENDIF
3150	!
3151	!
3152	! Compute the area of the basin
3153	!
3154	basin_area(ib,ij) = zero
3155	basin_hierarchy(ib,ij) = zero
3156	!
3157	SELECT CASE (hierar_method)
3158	!
3159	CASE("MINI")
3160	basin_hierarchy(ib,ij) = undef_sechiba
3161	!
3162	END SELECT
3163	basin_topoind(ib,ij) = zero
3164	!
3165	DO iz=1,basin_sz(ij)
3166	basin_area(ib,ij) = basin_area(ib,ij) + area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
3167	basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
3168	!
3169	! There are a number of ways to determine the hierarchy of the entire basin.
3170	! We allow for three here :
3171	! - Take the mean value
3172	! - Take the minimum value within the basin
3173	! - Take the value at the outflow point
3174	! Probably taking the value of the outflow point is the best solution.
3175	!
3176	SELECT CASE (hierar_method)
3177	!
3178	CASE("MEAN")
3179	! Mean hierarchy of the basin
3180	basin_hierarchy(ib,ij) = basin_hierarchy(ib,ij) + &
3181	& hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
3182	CASE("MINI")
3183	! The smalest value of the basin
3184	IF ( hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .LT. basin_hierarchy(ib,ij)) THEN
3185	basin_hierarchy(ib,ij) = hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
3186	ENDIF
3187	CASE("OUTP")
3188	! Value at the outflow point
3189	IF ( trip_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .GT. 100 ) THEN
3190	basin_hierarchy(ib,ij) = hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
3191	ENDIF
3192	CASE DEFAULT
3193	WRITE(numout,*) 'Unknown method for computing the hierarchy of the basin'
3194	STOP 'routing_globalize'
3195	END SELECT
3196	!
3197	ENDDO
3198	!
3199	basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std)
3200	!
3201	SELECT CASE (hierar_method)
3202	!
3203	CASE("MEAN")
3204	basin_hierarchy(ib,ij) = basin_hierarchy(ib,ij)/REAL(basin_sz(ij),r_std)
3205	!
3206	END SELECT
3207	!
3208	! To make sure that it has the lowest number if this is an outflow point we reset basin_hierarchy
3209	!
3210	IF (basin_bxout(ij) .LT. 0) THEN
3211	basin_hierarchy(ib,ij) = min_topoind
3212	basin_topoind(ib,ij) = min_topoind
3213	ENDIF
3214	!
3215	!
3216	! Keep the outflow boxes and basin
3217	!
3218	basin_flowdir(ib,ij) = basin_bxout(ij)
3219	IF (basin_bxout(ij) .GT. 0) THEN
3220	outflow_grid(ib,ij) = neighbours(ib,basin_bxout(ij))
3221	ELSE
3222	outflow_grid(ib,ij) = basin_bxout(ij)
3223	ENDIF
3224	!
3225	!
3226	ENDDO
3227	!
3228
3229	!
3230	END SUBROUTINE routing_globalize
3231	!
3232	!-----------------------------------------------------------------------------
3233	!
3234	SUBROUTINE routing_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, &
3235	& basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, nbcoastal,&
3236	& coastal_basin, invented_basins)
3237	!
3238	IMPLICIT NONE
3239	!
3240	! This subroutine will make the connections between the basins and ensure global coherence.
3241	!
3242	! The convention for outflow_grid is :
3243	! outflow_grid = -1 : River flow
3244	! outflow_grid = -2 : Coastal flow
3245	! outflow_grid = -3 : Return flow
3246	!
3247	!
3248	! INPUT
3249	!
3250	INTEGER(i_std), INTENT (in) :: nbpt
3251	INTEGER(i_std), DIMENSION(nbpt,8), INTENT (in) :: neighbours
3252	!
3253	INTEGER(i_std) :: nwbas
3254	INTEGER(i_std), DIMENSION(nbpt) :: basin_count
3255	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id
3256	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir
3257	REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area, basin_hierarchy
3258	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid
3259	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin
3260	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number
3261	INTEGER(i_std), DIMENSION(nbpt,nwbas,nbvmax) :: inflow_basin
3262	INTEGER(i_std), DIMENSION(nbpt,nwbas,nbvmax) :: inflow_grid
3263	INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal
3264	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin
3265	REAL(r_std), INTENT(in) :: invented_basins
3266	!
3267	! LOCAL
3268	!
3269	INTEGER(i_std) :: sp, sb, sbl, inp, bid, outdm1, outdp1
3270	INTEGER(i_std) :: dp1, dm1, dm1i, dp1i, bp1, bm1
3271	INTEGER(i_std) :: dop, bop
3272	INTEGER(i_std) :: fbas(nwbas), nbfbas
3273	REAL(r_std) :: fbas_hierarchy(nwbas)
3274	INTEGER(i_std) :: ff(1)
3275	!
3276	! WARNING
3277	LOGICAL, PARAMETER :: check = .FALSE.
3278	! ERRORS
3279	LOGICAL :: error1, error2, error3, error4, error5
3280
3281	error1=.FALSE.
3282	error2=.FALSE.
3283	error3=.FALSE.
3284	error4=.FALSE.
3285	error5=.FALSE.
3286
3287	outflow_basin(:,:) = undef_int
3288	inflow_number(:,:) = 0
3289	!
3290	DO sp=1,nbpt
3291	DO sb=1,basin_count(sp)
3292	!
3293	inp = outflow_grid(sp,sb)
3294	bid = basin_id(sp,sb)
3295	!
3296	! We only work on this point if it does not flow into the ocean
3297	! At this point any of the outflows is designated by a negative values in
3298	! outflow_grid
3299	!
3300	IF ( inp .GT. 0 ) THEN
3301	!
3302	! Now find the basin in the onflow point (inp)
3303	!
3304	nbfbas = 0
3305	!
3306	!
3307	DO sbl=1,basin_count(inp)
3308	!
3309	! Either it is a standard basin or one aggregated from ocean flow points.
3310	! If we flow into a another grid box we have to make sure that its hierarchy in the
3311	! basin is lower.
3312	!
3313	!
3314	IF ( basin_id(inp,sbl) .GT. 0 ) THEN
3315	IF ( basin_id(inp,sbl) .EQ. bid .OR. basin_id(inp,sbl) .GT. invented_basins) THEN
3316	nbfbas =nbfbas + 1
3317	fbas(nbfbas) = sbl
3318	fbas_hierarchy(nbfbas) = basin_hierarchy(inp,sbl)
3319	ENDIF
3320	ELSE
3321	IF ( COUNT(coastal_basin(inp,1:nbcoastal(inp)) .EQ. bid) .GT. 0 ) THEN
3322	nbfbas =nbfbas + 1
3323	fbas(nbfbas) = sbl
3324	fbas_hierarchy(nbfbas) = basin_hierarchy(inp,sbl)
3325	ENDIF
3326	ENDIF
3327	!
3328	ENDDO
3329	!
3330	! If we have more than one basin we will take the one which is lowest
3331	! in the hierarchy.
3332	!
3333	IF (nbfbas .GE. 1) THEN
3334	ff = MINLOC(fbas_hierarchy(1:nbfbas))
3335	sbl = fbas(ff(1))
3336	!
3337	bop = undef_int
3338	IF ( basin_hierarchy(inp,sbl) .LE. basin_hierarchy(sp,sb) ) THEN
3339	IF ( basin_hierarchy(inp,sbl) .LE. basin_hierarchy(sp,sb) ) THEN
3340	bop = sbl
3341	ELSE
3342	! The same hierarchy is allowed if both grids flow in about
3343	! the same direction :
3344	IF ( ( MOD(basin_flowdir(inp,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. &
3345	& ( basin_flowdir(inp,sbl) .EQ. basin_flowdir(sp,sb)).OR. &
3346	& ( MOD(basin_flowdir(inp,sbl)+7-1, 8)+1 .EQ. basin_flowdir(sp,sb)) ) THEN
3347	bop = sbl
3348	ENDIF
3349	ENDIF
3350	ENDIF
3351	!
3352	! If the basin is suitable (bop < undef_int) then take it
3353	!
3354	IF ( bop .LT. undef_int ) THEN
3355	outflow_basin(sp,sb) = bop
3356	inflow_number(inp,bop) = inflow_number(inp,bop) + 1
3357	IF ( inflow_number(inp,bop) .LE. nbvmax ) THEN
3358	inflow_grid(inp, bop, inflow_number(inp,bop)) = sp
3359	inflow_basin(inp, bop, inflow_number(inp,bop)) = sb
3360	ELSE
3361	error1=.TRUE.
3362	EXIT
3363	ENDIF
3364	ENDIF
3365	ENDIF
3366	!
3367	!
3368	ENDIF
3369	!
3370	!
3371	!
3372	! Did we find it ?
3373	!
3374	! In case the outflow point was ocean or we did not find the correct basin we start to look
3375	! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding
3376	! basin index (bp1 & bm1).
3377	!
3378	!
3379	IF ( outflow_basin(sp,sb) .EQ. undef_int &
3380	& .AND. basin_flowdir(sp,sb) .GT. 0) THEN
3381	!
3382	dp1i = MOD(basin_flowdir(sp,sb)+1-1, 8)+1
3383	dp1 = neighbours(sp,dp1i)
3384	dm1i = MOD(basin_flowdir(sp,sb)+7-1, 8)+1
3385	IF ( dm1i .LT. 1 ) dm1i = 8
3386	dm1 = neighbours(sp,dm1i)
3387	!
3388	!
3389	bp1 = -1
3390	IF ( dp1 .GT. 0 ) THEN
3391	DO sbl=1,basin_count(dp1)
3392	IF (basin_id(dp1,sbl) .EQ. bid .AND.&
3393	& basin_hierarchy(sp,sb) .GE. basin_hierarchy(dp1,sbl) .AND. &
3394	& bp1 .LT. 0) THEN
3395	IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dp1,sbl) ) THEN
3396	bp1 = sbl
3397	ELSE
3398	! The same hierarchy is allowed if both grids flow in about
3399	! the same direction :
3400	IF ( ( MOD(basin_flowdir(dp1,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. &
3401	& ( basin_flowdir(dp1,sbl) .EQ. basin_flowdir(sp,sb)).OR. &
3402	& ( MOD(basin_flowdir(dp1,sbl)+7-1, 8)+1 .EQ. basin_flowdir(sp,sb)) ) THEN
3403	bp1 = sbl
3404	ENDIF
3405	ENDIF
3406	ENDIF
3407	ENDDO
3408	ENDIF
3409	!
3410	bm1 = -1
3411	IF ( dm1 .GT. 0 ) THEN
3412	DO sbl=1,basin_count(dm1)
3413	IF (basin_id(dm1,sbl) .EQ. bid .AND.&
3414	& basin_hierarchy(sp,sb) .GE. basin_hierarchy(dm1,sbl) .AND. &
3415	& bm1 .LT. 0) THEN
3416	IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dm1,sbl) ) THEN
3417	bm1 = sbl
3418	ELSE
3419	! The same hierarchy is allowed if both grids flow in about
3420	! the same direction :
3421	IF ( ( MOD(basin_flowdir(dm1,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)) .OR. &
3422	& ( basin_flowdir(dm1,sbl) .EQ. basin_flowdir(sp,sb)) .OR. &
3423	& ( MOD(basin_flowdir(dm1,sbl)+7-1, 8)+1 .EQ. basin_flowdir(sp,sb)) ) THEN
3424	bm1 = sbl
3425	ENDIF
3426	ENDIF
3427	ENDIF
3428	ENDDO
3429	ENDIF
3430	!
3431	!
3432	! First deal with the case on land.
3433	!
3434	! For that we need to check if the water will be able to flow out of the grid dp1 or dm1
3435	! and not return to our current grid. If it is the current grid
3436	! then we can not do anything with that neighbour. Thus we set the
3437	! value of outdm1 and outdp1 back to -1
3438	!
3439	outdp1 = undef_int
3440	IF ( dp1 .GT. 0 .AND. bp1 .GT. 0 ) THEN
3441	! if the outflow is into the ocean then we put something less than undef_int in outdp1!
3442	IF (basin_flowdir(dp1,bp1) .GT. 0) THEN
3443	outdp1 = neighbours(dp1,basin_flowdir(dp1,bp1))
3444	IF ( outdp1 .EQ. sp ) outdp1 = undef_int
3445	ELSE
3446	outdp1 = nbpt + 1
3447	ENDIF
3448	ENDIF
3449	outdm1 = undef_int
3450	IF ( dm1 .GT. 0 .AND. bm1 .GT. 0 ) THEN
3451	IF (basin_flowdir(dm1,bm1) .GT. 0) THEN
3452	outdm1 = neighbours(dm1,basin_flowdir(dm1,bm1))
3453	IF ( outdm1 .EQ. sp ) outdm1 = undef_int
3454	ELSE
3455	outdm1 = nbpt + 1
3456	ENDIF
3457	ENDIF
3458	!
3459	! Now that we know our options we need go through them.
3460	!
3461	dop = undef_int
3462	bop = undef_int
3463	IF ( outdp1 .LT. undef_int .AND. outdm1 .LT. undef_int) THEN
3464	!
3465	! In this case we let the current basin flow into the smaller one
3466	!
3467	IF ( basin_area(dp1,bp1) .LT. basin_area(dm1,bm1) ) THEN
3468	dop = dp1
3469	bop = bp1
3470	ELSE
3471	dop = dm1
3472	bop = bm1
3473	ENDIF
3474	!
3475	!
3476	ELSE IF ( outdp1 .LT. undef_int ) THEN
3477	! If only the first one is possible
3478	dop = dp1
3479	bop = bp1
3480	ELSE IF ( outdm1 .LT. undef_int ) THEN
3481	! If only the second one is possible
3482	dop = dm1
3483	bop = bm1
3484	ELSE
3485	!
3486	! Now we are at the point where none of the neighboring points is suitable
3487	! or we have a coastal point.
3488	!
3489	! If there is an option to put the water into the ocean go for it.
3490	!
3491	IF ( outflow_grid(sp,sb) .LT. 0 .OR. dm1 .LT. 0 .OR. dp1 .LT. 0 ) THEN
3492	dop = -1
3493	ELSE
3494	!
3495	! If we are on a land point with only land neighbors but no one suitable to let the
3496	! water flow into we have to look for a solution in the current grid box.
3497	!
3498	!
3499	IF ( bp1 .LT. 0 .AND. bm1 .LT. 0 ) THEN
3500	!
3501	! Do we have more than one basin with the same ID ?
3502	!
3503	IF ( COUNT(basin_id(sp,1:basin_count(sp)) .EQ. bid) .GE. 2) THEN
3504	!
3505	! Now we can try the option of flowing into the basin of the same grid box.
3506	!
3507	DO sbl=1,basin_count(sp)
3508	IF (sbl .NE. sb .AND. basin_id(sp,sbl) .EQ. bid) THEN
3509	! In case this basin has a lower hierarchy or flows into a totaly
3510	! different direction we go for it.
3511	IF ( (basin_hierarchy(sp,sb) .GE. basin_hierarchy(sp,sbl)) .OR. &
3512	& (basin_flowdir(sp,sbl) .LT. dm1i .AND.&
3513	& basin_flowdir(sp,sbl) .GT. dp1i) ) THEN
3514	dop = sp
3515	bop = sbl
3516	IF (check) THEN
3517	IF (basin_hierarchy(sp,sb) .LT. basin_hierarchy(sp,sbl)) THEN
3518	WRITE(numout,*) '>>>>>>> POINT CORRECTED against hierarchy :',&
3519	& sp, sb, 'into', sbl
3520	ENDIF
3521	ENDIF
3522	ENDIF
3523	!
3524	ENDIF
3525	ENDDO
3526	!
3527	ENDIF
3528	ENDIF
3529	ENDIF
3530	!
3531	IF ( dop .EQ. undef_int .AND. bop .EQ. undef_int ) THEN
3532	IF (check) THEN
3533	WRITE(numout,*) 'Why are we here with point ', sp, sb
3534	WRITE(numout,*) 'Coodinates : (lon,lat) = ', lalo(sp,2), lalo(sp,1)
3535	WRITE(numout,*) 'Contfrac : = ', contfrac(sp)
3536	WRITE(numout,*) 'Local Basin ID :', basin_id(sp,1:basin_count(sp))
3537	WRITE(numout,*) 'Local hierarchies :', basin_hierarchy(sp,1:basin_count(sp))
3538	WRITE(numout,*) 'Local basin_flowdir :', basin_flowdir(sp,1:basin_count(sp))
3539	WRITE(numout,*) 'Local outflowgrid :', outflow_grid(sp,1:basin_count(sp))
3540	WRITE(numout,*) 'outflow_grid :', inp
3541	WRITE(numout,*) 'Coodinates outflow : (lon,lat) = ', lalo(inp,2), lalo(inp,1)
3542	WRITE(numout,*) 'Contfrac : = ', contfrac(inp)
3543	WRITE(numout,*) 'Outflow Basin ID :', basin_id(inp,1:basin_count(inp))
3544	WRITE(numout,*) 'Outflow hierarchies :', basin_hierarchy(inp,1:basin_count(inp))
3545	WRITE(numout,*) 'Outflow basin_flowdir :', basin_flowdir(inp,1:basin_count(inp))
3546	WRITE(numout,*) 'Explored options +1 :', dp1, bp1, outdp1
3547	WRITE(numout,*) 'Explored +1 Basin ID :', basin_id(dp1,1:basin_count(dp1))
3548	WRITE(numout,*) 'Explored +1 hierarchies :', basin_hierarchy(dp1,1:basin_count(dp1))
3549	WRITE(numout,*) 'Explored +1 basin_flowdir :', basin_flowdir(dp1,1:basin_count(dp1))
3550	WRITE(numout,*) 'Explored options -1 :', dm1, bm1, outdm1
3551	WRITE(numout,*) 'Explored -1 Basin ID :', basin_id(dm1,1:basin_count(dm1))
3552	WRITE(numout,*) 'Explored -1 hierarchies :', basin_hierarchy(dm1,1:basin_count(dm1))
3553	WRITE(numout,*) 'Explored -1 basin_flowdir :', basin_flowdir(dm1,1:basin_count(dm1))
3554	WRITE(numout,) '***************************'
3555	ENDIF
3556	IF ( contfrac(sp) > 0.01 ) THEN
3557	error2=.TRUE.
3558	EXIT
3559	ENDIF
3560	ENDIF
3561	!
3562	ENDIF
3563	!
3564	! Now that we know where we want the water to flow to we write the
3565	! the information in the right fields.
3566	!
3567	IF ( dop .GT. 0 .AND. dop .NE. undef_int ) THEN
3568	outflow_grid(sp,sb) = dop
3569	outflow_basin(sp,sb) = bop
3570	inflow_number(dop,bop) = inflow_number(dop,bop) + 1
3571	IF ( inflow_number(dop,bop) .LE. nbvmax ) THEN
3572	inflow_grid(dop, bop, inflow_number(dop,bop)) = sp
3573	inflow_basin(dop, bop, inflow_number(dop,bop)) = sb
3574	ELSE
3575	error3=.TRUE.
3576	EXIT
3577	ENDIF
3578	!
3579	ELSE
3580	outflow_grid(sp,sb) = -2
3581	outflow_basin(sp,sb) = undef_int
3582	ENDIF
3583	!
3584	ENDIF
3585	!
3586	!
3587	! If we still have not found anything then we have to check that there is not a basin
3588	! within the same grid box which has a lower hierarchy.
3589	!
3590	!
3591	IF ( outflow_grid(sp,sb) .GT. 0 .AND. outflow_basin(sp,sb) .EQ. undef_int &
3592	& .AND. basin_flowdir(sp,sb) .GT. 0) THEN
3593	!
3594
3595	IF (check) &
3596	WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb
3597	!
3598	DO sbl=1,basin_count(sp)
3599	!
3600	! Three conditons are needed to let the water flow into another basin of the
3601	! same grid :
3602	! - another basin than the current one
3603	! - same ID
3604	! - of lower hierarchy.
3605	!
3606	IF ( (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid)&
3607	& .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl)) ) THEN
3608	outflow_basin(sp,sb) = sbl
3609	inflow_number(sp,sbl) = inflow_number(sp,sbl) + 1
3610	IF ( inflow_number(sp,sbl) .LE. nbvmax ) THEN
3611	IF ( sp .EQ. 42 .AND. sbl .EQ. 1) THEN
3612	IF (check) &
3613	WRITE(numout,*) 'ADD INFLOW (3):', sp, sb
3614	ENDIF
3615	inflow_grid(sp, sbl, inflow_number(sp,sbl)) = sp
3616	inflow_basin(sp, sbl, inflow_number(sp,sbl)) = sb
3617	ELSE
3618	error4=.TRUE.
3619	EXIT
3620	ENDIF
3621	ENDIF
3622	ENDDO
3623	ENDIF
3624	!
3625	! Ok that is it, we give up :-)
3626	!
3627	IF ( outflow_grid(sp,sb) .GT. 0 .AND. outflow_basin(sp,sb) .EQ. undef_int &
3628	& .AND. basin_flowdir(sp,sb) .GT. 0) THEN
3629	!
3630	error5=.TRUE.
3631	EXIT
3632	ENDIF
3633	ENDDO
3634	!
3635	ENDDO
3636	IF (error1) THEN
3637	WRITE(numout,*) " routing_linkup : bop .LT. undef_int",bop
3638	CALL ipslerr(3,'routing_linkup', &
3639	"bop .LT. undef_int",'Increase nbvmax','stop routing_linkup')
3640	ENDIF
3641	IF (error2) THEN
3642	CALL ipslerr(3,'routing_linkup', &
3643	& 'In the routine which make connections between the basins and ensure global coherence,', &
3644	& 'there is a problem with outflow linkup without any valid direction. Try with check=.TRUE.', &
3645	& '(Perhaps there is a problem with the grid.)')
3646	ENDIF
3647	IF (error3) THEN
3648	WRITE(numout,*) " routing_linkup : dop .GT. 0 .AND. dop .NE. undef_int",dop
3649	CALL ipslerr(3,'routing_linkup', &
3650	"dop .GT. 0 .AND. dop .NE. undef_int",'Increase nbvmax. Try with check=.TRUE.','stop routing_linkup')
3651	ENDIF
3652	IF (error4) THEN
3653	WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", &
3654	& " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, &
3655	& basin_hierarchy(sp,sb),basin_hierarchy(sp,sbl)
3656	CALL ipslerr(3,'routing_linkup', &
3657	"(sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))" &
3658	,'Increase nbvmax. Try with check=.TRUE.','stop routing_linkup')
3659	ENDIF
3660	IF (error5) THEN
3661	WRITE(numout,*) 'We could not find the basin into which we need to flow'
3662	WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb
3663	WRITE(numout,*) 'Explored neighbours :', dm1, dp1
3664	WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb)
3665	WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb)
3666	WRITE(numout,*) 'Outlfow basin :',outflow_basin(sp,sb)
3667	WRITE(numout,*) 'basin ID:',basin_id(sp,sb)
3668	WRITE(numout,*) 'Hierarchy :', basin_hierarchy(sp,sb)
3669	CALL ipslerr(3,'routing_linkup', &
3670	"We could not find the basin into which we need to flow",'Try with check=.TRUE.','stop routing_linkup')
3671	ENDIF
3672	!
3673	! Check for each outflow basin that it exists
3674	!
3675	DO sp=1,nbpt
3676	DO sb=1,basin_count(sp)
3677	!
3678	inp = outflow_grid(sp,sb)
3679	sbl = outflow_basin(sp,sb)
3680	IF ( inp .GE. 0 ) THEN
3681	IF ( basin_count(inp) .LT. sbl ) THEN
3682	WRITE(numout,*) 'point :', sp, ' basin :', sb
3683	WRITE(numout,*) 'Flows into point :', inp, ' basin :', sbl
3684	WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(inp)
3685	STOP
3686	ENDIF
3687	ENDIF
3688	ENDDO
3689	ENDDO
3690	!
3691	END SUBROUTINE routing_linkup
3692	!
3693	!---------------------------------------------------------------------------
3694	!
3695	SUBROUTINE routing_fetch(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_id,&
3696	& outflow_grid, outflow_basin, fetch_basin)
3697	!
3698	IMPLICIT NONE
3699	!
3700	! This subroutine will compute the fetch of each basin. This means that for each basin we
3701	! will know how much area is upstream. It will help decide how to procede with the
3702	! the truncation later and allow to set correctly in outflow_grid the distinction
3703	! between coastal and river flow
3704	!
3705	!
3706	! INPUT
3707	!
3708	INTEGER(i_std), INTENT(in) :: nbpt
3709	!
3710	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution
3711	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of land in each grid box
3712	!
3713	INTEGER(i_std) :: nwbas
3714	INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count
3715	REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area
3716	INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id
3717	INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid
3718	INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin
3719	REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: fetch_basin
3720	!
3721	! LOCAL
3722	!
3723	INTEGER(i_std) :: ib, ij, ff(1), it, itt, igrif, ibasf, nboutflow
3724	REAL(r_std) :: contarea, totbasins
3725	REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area
3726	INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex
3727	!
3728	!
3729	! Normalize the area of all basins
3730	!
3731	DO ib=1,nbpt
3732	!
3733	totbasins = SUM(basin_area(ib,1:basin_count(ib)))
3734	contarea = resolution(ib,1)resolution(ib,2)contfrac(ib)
3735	!
3736	DO ij=1,basin_count(ib)
3737	basin_area(ib,ij) = basin_area(ib,ij)/totbasins*contarea
3738	ENDDO
3739	!
3740	ENDDO
3741	WRITE(numout,*) 'Normalization done'
3742	!
3743	! Compute the area upstream of each basin
3744	!
3745	fetch_basin(:,:) = zero
3746	!
3747	!
3748	DO ib=1,nbpt
3749	!
3750	DO ij=1,basin_count(ib)
3751	!
3752	fetch_basin(ib, ij) = fetch_basin(ib, ij) + basin_area(ib,ij)
3753	!
3754	igrif = outflow_grid(ib,ij)
3755	ibasf = outflow_basin(ib,ij)
3756	!
3757	itt = 0
3758	DO WHILE (igrif .GT. 0)
3759	fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ib, ij)
3760	it = outflow_grid(igrif, ibasf)
3761	ibasf = outflow_basin(igrif, ibasf)
3762	igrif = it
3763	itt = itt + 1
3764	IF ( itt .GT. 500) THEN
3765	WRITE(numout,&
3766	"('Grid ',I5, ' and basin ',I5, 'did not converge after iteration ',I5)") ib, ij, itt
3767	WRITE(numout,*) 'Basin ID :', basin_id(igrif,ibasf)
3768	WRITE(numout,&
3769	"('We are stuck with the flow into grid ',I5,' and basin ',I5)") igrif, ibasf
3770	IF ( itt .GT. 510) THEN
3771	STOP ' routing_fetch'
3772	ENDIF
3773	ENDIF
3774	ENDDO
3775	!
3776	ENDDO
3777	!
3778	ENDDO
3779	!
3780	WRITE(numout,*) 'The smalest FETCH :', MINVAL(fetch_basin)
3781	WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin)
3782	!
3783	! Now we set for the 'num_largest' largest basins the outflow condition as stream flow
3784	! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow
3785	! (i.e. outflow_grid = -2). The return flow is not touched.
3786	!
3787	nboutflow = 0
3788	!
3789	DO ib=1,nbpt
3790	!
3791	DO ij=1,basin_count(ib)
3792	!
3793	! We do not need any more the river flow flag as we are going to reset it.
3794	!
3795	IF ( outflow_grid(ib,ij) .EQ. -1) THEN
3796	outflow_grid(ib,ij) = -2
3797	ENDIF
3798	!
3799	IF ( outflow_grid(ib,ij) .EQ. -2) THEN
3800	!
3801	nboutflow = nboutflow + 1
3802	tmp_area(nboutflow) = fetch_basin(ib,ij)
3803	tmpindex(nboutflow,1) = ib
3804	tmpindex(nboutflow,2) = ij
3805	!
3806	ENDIF
3807	!
3808	ENDDO
3809	ENDDO
3810	!
3811	DO ib=1, num_largest
3812	ff = MAXLOC(tmp_area(1:nboutflow))
3813	outflow_grid(tmpindex(ff(1),1), tmpindex(ff(1),2)) = -1
3814	tmp_area(ff(1)) = zero
3815	ENDDO
3816	!
3817	END SUBROUTINE routing_fetch
3818	!
3819	!---------------------------------------------------------------------------
3820	!
3821	SUBROUTINE routing_truncate(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_topoind,&
3822	& fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
3823	& inflow_grid, inflow_basin)
3824	!
3825	IMPLICIT NONE
3826	!
3827	!
3828	! This subroutine reduces the number of basins per gird to the value chosen by the user.
3829	! it also computes the final field which will be used to route the water at the
3830	! requested truncation.
3831	!
3832	! INPUT
3833	!
3834	INTEGER(i_std) :: nbpt
3835	!
3836	REAL(r_std), DIMENSION(nbpt,2) :: resolution
3837	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of land in each grid box
3838	!
3839	INTEGER(i_std) :: nwbas
3840	INTEGER(i_std), DIMENSION(nbpt) :: basin_count
3841	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id
3842	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir
3843	REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area
3844	REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind
3845	REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin
3846	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid
3847	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin
3848	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number
3849	INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_basin
3850	INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_grid
3851	!
3852	! LOCAL
3853	!
3854	INTEGER(i_std), PARAMETER :: pickmax = 200
3855	INTEGER(i_std) :: ib, ij, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2)
3856	INTEGER(i_std) :: ii, kbas, sbas, ik, iter, ibt, obj
3857	REAL(r_std), DIMENSION(nbpt,nbasmax) :: floflo
3858	REAL(r_std), DIMENSION(nbpt) :: gridarea, gridbasinarea
3859	REAL(r_std) :: ratio
3860	INTEGER(i_std), DIMENSION(pickmax,2) :: largest_basins
3861	INTEGER(i_std), DIMENSION(pickmax) :: tmp_ids
3862	INTEGER(i_std) :: multbas, iml(1)
3863	INTEGER(i_std), DIMENSION(pickmax) :: multbas_sz
3864	REAL(r_std), DIMENSION(pickmax) :: tmp_area
3865	INTEGER(i_std), DIMENSION(pickmax,pickmax) :: multbas_list
3866	!
3867	INTEGER(i_std) :: nbtruncate
3868	INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: indextrunc
3869	!
3870	IF ( .NOT. ALLOCATED(indextrunc)) THEN
3871	ALLOCATE(indextrunc(nbpt))
3872	ENDIF
3873	!
3874	! Truncate if needed and find the path closest to the high resolution data.
3875	!
3876	! The algorithm :
3877	!
3878	! We only go through this procedure only as many times as there are basins to take out at most.
3879	! This is important as it allows the simplifications to spread from one grid to the other.
3880	! The for each step of the iteration and at each grid point we check the following options for
3881	! simplifying the pathways of water :
3882	! 1) If the basin of a grid flows into another basin of the same grid. Kill the one which only
3883	! served as a transition
3884	! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow.
3885	! We kill the smallest one and put into the largest basin. There is no need to manage many
3886	! basins going into the ocean as coastal flows.
3887	! 3) If we have streams run in parallel from one gird box to the others (that is these are
3888	! different basins) we will put the smaller one in the larger one. This may hapen at any
3889	! level of the flow but in theory it should propagate downstream.
3890	! 4) If we have two basins with the same ID but flow into different grid boxes we sacrifice
3891	! the smallest one and route it through the largest.
3892	!
3893	! Obviously if any of the options find something then we skip the rest and take out the basin.
3894	!
3895	!
3896	! We have to go through the grid as least as often as we have to reduce the number of basins
3897	! For good measure we add 3 more passages.
3898	!
3899	!
3900	DO iter = 1, MAXVAL(basin_count) - nbasmax +3
3901	!
3902	! Get the points over which we wish to truncate
3903	!
3904	nbtruncate = 0
3905	DO ib = 1, nbpt
3906	IF ( basin_count(ib) .GT. nbasmax ) THEN
3907	nbtruncate = nbtruncate + 1
3908	indextrunc(nbtruncate) = ib
3909	ENDIF
3910	ENDDO
3911	!
3912	! Go through the basins which need to be truncated.
3913	!
3914	DO ibt=1,nbtruncate
3915	!
3916	ib = indextrunc(ibt)
3917	!
3918	! Check if we have basin which flows into a basin in the same grid
3919	! kbas = basin we will have to kill
3920	! sbas = basin which takes over kbas
3921	!
3922	kbas = 0
3923	sbas = 0
3924	!
3925	! 1) Can we find a basin which flows into a basin of the same grid ?
3926	!
3927	DO ij=1,basin_count(ib)
3928	DO ii=1,basin_count(ib)
3929	IF ( outflow_grid(ib,ii) .EQ. ib .AND. outflow_basin(ib, ii) .EQ. ij .AND. kbas*sbas .NE. 0) THEN
3930	kbas = ii
3931	sbas = ij
3932	ENDIF
3933	ENDDO
3934	ENDDO
3935	!
3936	! 2) Merge two basins which flow into the ocean as coastal or return flow
3937	! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and
3938	! have not found anything yet!
3939	!
3940	IF ( (COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 .OR.&
3941	& COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -3) .GT. 1) .AND.&
3942	& kbas*sbas .EQ. 0) THEN
3943	!
3944	multbas = 0
3945	multbas_sz(:) = 0
3946	!
3947	IF ( COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 ) THEN
3948	obj = -2
3949	ELSE
3950	obj = -3
3951	ENDIF
3952	!
3953	! First we get the list of all basins which go out as coastal or return flow (obj)
3954	!
3955	DO ij=1,basin_count(ib)
3956	IF ( outflow_grid(ib,ij) .EQ. obj ) THEN
3957	multbas = multbas + 1
3958	multbas_sz(multbas) = ij
3959	tmp_area(multbas) = fetch_basin(ib,ij)
3960	ENDIF
3961	ENDDO
3962	!
3963	! Now the take the smalest to be transfered to the largest
3964	!
3965	iml = MAXLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
3966	sbas = multbas_sz(iml(1))
3967	iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
3968	kbas = multbas_sz(iml(1))
3969	!
3970	ENDIF
3971	!
3972	! 3) If we have basins flowing into the same grid but different basins then we put them
3973	! together. Obviously we first work with the grid which has most streams runing into it
3974	! and puting the smalest in the largests catchments.
3975	!
3976	IF ( kbas*sbas .EQ. 0) THEN
3977	!
3978	tmp_ids(1:basin_count(ib)) = outflow_grid(ib,1:basin_count(ib))
3979	multbas = 0
3980	multbas_sz(:) = 0
3981	!
3982	! First obtain the list of basins which flow into the same basin
3983	!
3984	DO ij=1,basin_count(ib)
3985	IF ( outflow_grid(ib,ij) .GT. 0 .AND.&
3986	& COUNT(tmp_ids(1:basin_count(ib)) .EQ. outflow_grid(ib,ij)) .GT. 1) THEN
3987	multbas = multbas + 1
3988	DO ii=1,basin_count(ib)
3989	IF ( tmp_ids(ii) .EQ. outflow_grid(ib,ij)) THEN
3990	multbas_sz(multbas) = multbas_sz(multbas) + 1
3991	multbas_list(multbas,multbas_sz(multbas)) = ii
3992	tmp_ids(ii) = -99
3993	ENDIF
3994	ENDDO
3995	ELSE
3996	tmp_ids(ij) = -99
3997	ENDIF
3998	ENDDO
3999	!
4000	! Did we come up with any basins to deal with this way ?
4001	!
4002	IF ( multbas .GT. 0 ) THEN
4003	!
4004	iml = MAXLOC(multbas_sz(1:multbas))
4005	ik = iml(1)
4006	!
4007	! Take the smalest and largest of these basins !
4008	!
4009	DO ii=1,multbas_sz(ik)
4010	tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
4011	ENDDO
4012	iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
4013	sbas = multbas_list(ik,iml(1))
4014	iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
4015	kbas = multbas_list(ik,iml(1))
4016	!
4017	ENDIF
4018	!
4019	ENDIF
4020	!
4021	! 4) If we have twice the same basin we put them together even if they flow into different
4022	! directions. If one of them goes to the ocean it takes the advantage.
4023	!
4024	IF ( kbas*sbas .EQ. 0) THEN
4025	!
4026	tmp_ids(1:basin_count(ib)) = basin_id(ib,1:basin_count(ib))
4027	multbas = 0
4028	multbas_sz(:) = 0
4029	!
4030	! First obtain the list of basins which have sub-basins in this grid box.
4031	! (these are identified by their IDs)
4032	!
4033	DO ij=1,basin_count(ib)
4034	IF ( COUNT(tmp_ids(1:basin_count(ib)) .EQ. basin_id(ib,ij)) .GT. 1) THEN
4035	multbas = multbas + 1
4036	DO ii=1,basin_count(ib)
4037	IF ( tmp_ids(ii) .EQ. basin_id(ib,ij)) THEN
4038	multbas_sz(multbas) = multbas_sz(multbas) + 1
4039	multbas_list(multbas,multbas_sz(multbas)) = ii
4040	tmp_ids(ii) = -99
4041	ENDIF
4042	ENDDO
4043	ELSE
4044	tmp_ids(ij) = -99
4045	ENDIF
4046	ENDDO
4047	!
4048	! We are going to work on the basin with the largest number of sub-basins.
4049	! (IF we have a basin which has subbasins !)
4050	!
4051	IF ( multbas .GT. 0 ) THEN
4052	!
4053	iml = MAXLOC(multbas_sz(1:multbas))
4054	ik = iml(1)
4055	!
4056	! If one of the basins goes to the ocean then it is going to have the priority
4057	!
4058	tmp_area(:) = zero
4059	IF ( COUNT(outflow_grid(ib,multbas_list(ik,1:multbas_sz(ik))) .LT. 0) .GT. 0) THEN
4060	DO ii=1,multbas_sz(ik)
4061	IF ( outflow_grid(ib,multbas_list(ik,ii)) .LT. 0 .AND. sbas .EQ. 0 ) THEN
4062	sbas = multbas_list(ik,ii)
4063	ELSE
4064	tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
4065	ENDIF
4066	ENDDO
4067	! take the smalest of the subbasins
4068	iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
4069	kbas = multbas_list(ik,iml(1))
4070	ELSE
4071	!
4072	! Else we take simply the largest and smalest
4073	!
4074	DO ii=1,multbas_sz(ik)
4075	tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
4076	ENDDO
4077	iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
4078	sbas = multbas_list(ik,iml(1))
4079	iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
4080	kbas = multbas_list(ik,iml(1))
4081	!
4082	ENDIF
4083	!
4084	!
4085	ENDIF
4086	ENDIF
4087	!
4088	!
4089	!
4090	! Then we call routing_killbas to clean up the basins in this grid
4091	!
4092	IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
4093	CALL routing_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,&
4094	& fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
4095	& inflow_grid, inflow_basin)
4096	ENDIF
4097	!
4098	ENDDO
4099	!
4100	!
4101	ENDDO
4102	!
4103	! If there are any grids left with too many basins we need to take out the big hammer !
4104	! We will only do it if this represents less than 5% of all points.
4105	!
4106	IF ( COUNT(basin_count .GT. nbasmax) .GT. 0 ) THEN
4107	!
4108	!
4109	IF ( COUNT(basin_count .GT. nbasmax)/nbpt*100 .GT. 5 ) THEN
4110	WRITE(numout,) 'We have ', COUNT(basin_count .GT. nbasmax)/nbpt100, '% of all points which do not yet'
4111	WRITE(numout,*) 'have the right trunctaction. That is too much to apply a brutal method'
4112	DO ib = 1, nbpt
4113	IF ( basin_count(ib) .GT. nbasmax ) THEN
4114	!
4115	WRITE(numout,*) 'We did not find a basin which could be supressed. We will'
4116	WRITE(numout,*) 'not be able to reduce the truncation in grid ', ib
4117	DO ij=1,basin_count(ib)
4118	WRITE(numout,*) 'grid, basin nb and id :', ib, ij, basin_id(ib,ij)
4119	WRITE(numout,*) 'Outflow grid and basin ->', outflow_grid(ib,ij), outflow_basin(ib, ij)
4120	ENDDO
4121	ENDIF
4122	ENDDO
4123	STOP 'routing_truncate'
4124	!
4125	ELSE
4126	!
4127	!
4128	DO ib = 1,nbpt
4129	DO WHILE ( basin_count(ib) .GT. nbasmax )
4130	!
4131	WRITE(numout,*) 'HAMMER, ib, basin_count :', ib, basin_count(ib)
4132	!
4133	! Here we simply put the smallest basins into the largest ones. It is really a brute force
4134	! method but it will only be applied if everything has failed.
4135	!
4136	DO ii = 1,basin_count(ib)
4137	tmp_area(ii) = fetch_basin(ib, ii)
4138	ENDDO
4139	!
4140	iml = MAXLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
4141	sbas =iml(1)
4142	iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
4143	kbas = iml(1)
4144	!
4145	IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
4146	CALL routing_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,&
4147	& fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
4148	& inflow_grid, inflow_basin)
4149	ENDIF
4150	ENDDO
4151	ENDDO
4152	!
4153	ENDIF
4154	!
4155	!
4156	ENDIF
4157	!
4158	! Now that we have reached the right truncation (resolution) we will start
4159	! to produce the variables we will use to route the water.
4160	!
4161	DO ib=1,nbpt
4162	!
4163	! For non existing basins the route_tobasin variable is put to zero. This will allow us
4164	! to pick up the number of basin afterwards.
4165	!
4166	route_togrid(ib,:) = ib
4167	route_tobasin(ib,:) = 0
4168	routing_area(ib,:) = zero
4169	!
4170	ENDDO
4171	!
4172	! Transfer the info into the definitive variables
4173	!
4174	DO ib=1,nbpt
4175	DO ij=1,basin_count(ib)
4176	routing_area(ib,ij) = basin_area(ib,ij)
4177	topo_resid(ib,ij) = basin_topoind(ib,ij)
4178	global_basinid(ib,ij) = basin_id(ib,ij)
4179	route_togrid(ib,ij) = outflow_grid(ib,ij)
4180	route_tobasin(ib,ij) = outflow_basin(ib,ij)
4181	ENDDO
4182	ENDDO
4183	!
4184	!
4185	! Set the new convention for the outflow conditions
4186	! Now it is based in the outflow basin and the outflow grid will
4187	! be the same as the current.
4188	! returnflow to the grid : nbasmax + 1
4189	! coastal flow : nbasmax + 2
4190	! river outflow : nbasmax + 3
4191	!
4192	! Here we put everything here in coastal flow. It is later where we will
4193	! put the largest basins into river outflow.
4194	!
4195	DO ib=1,nbpt
4196	DO ij=1,basin_count(ib)
4197	! River flows
4198	IF ( route_togrid(ib,ij) .EQ. -1 ) THEN
4199	route_tobasin(ib,ij) = nbasmax + 2
4200	route_togrid(ib,ij) = ib
4201	! Coastal flows
4202	ELSE IF ( route_togrid(ib,ij) .EQ. -2 ) THEN
4203	route_tobasin(ib,ij) = nbasmax + 2
4204	route_togrid(ib,ij) = ib
4205	! Return flow
4206	ELSE IF ( route_togrid(ib,ij) .EQ. -3 ) THEN
4207	route_tobasin(ib,ij) = nbasmax + 1
4208	route_togrid(ib,ij) = ib
4209	ENDIF
4210	ENDDO
4211	ENDDO
4212	!
4213	! A second check on the data. Just make sure that each basin flows somewhere.
4214	!
4215	DO ib=1,nbpt
4216	DO ij=1,basin_count(ib)
4217	ibf = route_togrid(ib,ij)
4218	ijf = route_tobasin(ib,ij)
4219	IF ( ijf .GT. basin_count(ibf) .AND. ijf .LE. nbasmax) THEN
4220	WRITE(numout,*) 'Second check'
4221	WRITE(numout,*) 'point :', ib, ' basin :', ij
4222	WRITE(numout,*) 'Flows into point :', ibf, ' basin :', ijf
4223	WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(ibf)
4224	STOP
4225	ENDIF
4226	ENDDO
4227	ENDDO
4228	!
4229	! Verify areas of the contienents
4230	!
4231	floflo(:,:) = zero
4232	gridarea(:) = contfrac(:)resolution(:,1)resolution(:,2)
4233	DO ib=1,nbpt
4234	gridbasinarea(ib) = SUM(routing_area(ib,:))
4235	ENDDO
4236	!
4237	DO ib=1,nbpt
4238	DO ij=1,basin_count(ib)
4239	cnt = 0
4240	igrif = ib
4241	ibasf = ij
4242	DO WHILE (ibasf .LE. nbasmax .AND. cnt .LT. nbasmax*nbpt)
4243	cnt = cnt + 1
4244	pold = igrif
4245	bold = ibasf
4246	igrif = route_togrid(pold, bold)
4247	ibasf = route_tobasin(pold, bold)
4248	IF ( ibasf .GT. basin_count(igrif) .AND. ibasf .LE. nbasmax) THEN
4249	WRITE(numout,*) 'We should not be here as the basin flows into the pampa'
4250	WRITE(numout,*) 'Last correct point :', pold, bold
4251	WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold)
4252	WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold)
4253	WRITE(numout,*) 'Where we ended up :', igrif,ibasf
4254	STOP
4255	ENDIF
4256	ENDDO
4257	!
4258	IF ( ibasf .GT. nbasmax ) THEN
4259	floflo(igrif,bold) = floflo(igrif,bold) + routing_area(ib,ij)
4260	ELSE
4261	WRITE(numout,*) 'The flow did not end up in the ocean or in the grid cell.'
4262	WRITE(numout,*) 'For grid ', ib, ' and basin ', ij
4263	WRITE(numout,*) 'The last grid was ', igrif, ' and basin ', ibasf
4264	STOP 'routing_truncate'
4265	ENDIF
4266	ENDDO
4267	ENDDO
4268	!
4269	DO ib=1,nbpt
4270	IF ( gridbasinarea(ib) > zero ) THEN
4271	ratio = gridarea(ib)/gridbasinarea(ib)
4272	routing_area(ib,:) = routing_area(ib,:)*ratio
4273	ENDIF
4274	ENDDO
4275	!
4276	WRITE(numout,*) 'Sum of area of all outflow areas :',SUM(routing_area)
4277	WRITE(numout,*) 'Surface of all continents :', SUM(gridarea)
4278	!
4279	! Redo the the distinction between river outflow and coastal flow. We can not
4280	! take into account the return flow points.
4281	!
4282	ibf = 0
4283	DO ib=1, pickmax
4284	ff = MAXLOC(floflo)
4285	IF ( route_tobasin(ff(1), ff(2)) .EQ. nbasmax + 2) THEN
4286	ibf = ibf + 1
4287	largest_basins(ibf,:) = ff(:)
4288	ENDIF
4289	floflo(ff(1), ff(2)) = zero
4290	ENDDO
4291	!
4292	! Put the largest basins into river flows.
4293	!
4294	IF ( ibf .LT. num_largest) THEN
4295	WRITE(numout,*) 'Not enough basins to choose the ', num_largest, 'largest'
4296	STOP 'routing_truncate'
4297	ENDIF
4298	!
4299	!
4300	!
4301	DO ib=1, num_largest
4302	route_tobasin(largest_basins(ib,1),largest_basins(ib,2)) = nbasmax + 3
4303	ENDDO
4304	!
4305	WRITE(numout,*) 'NUMBER OF RIVERS :', COUNT(route_tobasin .GE. nbasmax + 3)
4306	!
4307	END SUBROUTINE routing_truncate
4308	!
4309	!-------------------------------------------------------------------------------------
4310	!
4311	SUBROUTINE routing_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, basin_topoind,&
4312	& fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
4313	& inflow_grid, inflow_basin)
4314	!
4315	!
4316	IMPLICIT NONE
4317	!
4318	! The aim of this routine is to kill a basin (that is put into another larger one). When
4319	! we do this we need to be careful and change all associated variables.
4320	!
4321	INTEGER(i_std) :: tokill, totakeover
4322	INTEGER(i_std) :: nbpt, ib
4323	!
4324	INTEGER(i_std) :: nwbas
4325	INTEGER(i_std), DIMENSION(nbpt) :: basin_count
4326	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id
4327	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir
4328	REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area
4329	REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind
4330	REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin
4331	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid
4332	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin
4333	INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number
4334	INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_basin
4335	INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_grid
4336	!
4337	! LOCAL
4338	!
4339	INTEGER(i_std) :: inf, ibs, ing, inb, ibasf, igrif, it
4340	LOGICAL :: doshift
4341	!
4342	! Update the information needed in the basin "totakeover"
4343	! For the moment only area
4344	!!$ !
4345	!!$ WRITE(numout,*) 'KILL BASIN :', ib, tokill, totakeover, basin_id(ib,tokill), basin_id(ib,totakeover)
4346	!!$ !
4347	!
4348	basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib, tokill)
4349	basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover) + basin_topoind(ib, tokill))/2.0
4350	!
4351	! Add the fetch of the basin will kill to the one which gets the water
4352	!
4353	fetch_basin(ib, totakeover) = fetch_basin(ib, totakeover) + fetch_basin(ib, tokill)
4354	igrif = outflow_grid(ib,totakeover)
4355	ibasf = outflow_basin(ib,totakeover)
4356	!
4357	inf = 0
4358	DO WHILE (igrif .GT. 0)
4359	fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill)
4360	it = outflow_grid(igrif, ibasf)
4361	ibasf = outflow_basin(igrif, ibasf)
4362	igrif = it
4363	inf = inf + 1
4364	ENDDO
4365	!
4366	! Take out the basin we have just rerouted from the fetch of the basins in which it used to flow.
4367	!
4368	igrif = outflow_grid(ib,tokill)
4369	ibasf = outflow_basin(ib,tokill)
4370	!
4371	DO WHILE (igrif .GT. 0)
4372	fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) - fetch_basin(ib, tokill)
4373	it = outflow_grid(igrif, ibasf)
4374	ibasf = outflow_basin(igrif, ibasf)
4375	igrif = it
4376	ENDDO
4377	!
4378	! Redirect the flows which went into the basin to be killed before we change everything
4379	!
4380	DO inf = 1, inflow_number(ib, tokill)
4381	outflow_basin(inflow_grid(ib, tokill, inf), inflow_basin(ib, tokill, inf)) = totakeover
4382	inflow_number(ib, totakeover) = inflow_number(ib, totakeover) + 1
4383	inflow_grid(ib, totakeover, inflow_number(ib, totakeover)) = inflow_grid(ib, tokill, inf)
4384	inflow_basin(ib, totakeover, inflow_number(ib, totakeover)) = inflow_basin(ib, tokill, inf)
4385	ENDDO
4386	!
4387	! Take out the basin to be killed from the list of inflow basins of the downstream basin
4388	! (In case the basin does not flow into an ocean or lake)
4389	!
4390	IF ( outflow_grid(ib,tokill) .GT. 0) THEN
4391	!
4392	ing = outflow_grid(ib, tokill)
4393	inb = outflow_basin(ib, tokill)
4394	doshift = .FALSE.
4395	!
4396	DO inf = 1, inflow_number(ing, inb)
4397	IF ( doshift ) THEN
4398	inflow_grid(ing, inb, inf-1) = inflow_grid(ing, inb, inf)
4399	inflow_basin(ing, inb, inf-1) = inflow_basin(ing, inb, inf)
4400	ENDIF
4401	IF ( inflow_grid(ing, inb, inf) .EQ. ib .AND. inflow_basin(ing, inb, inf) .EQ. tokill) THEN
4402	doshift = .TRUE.
4403	ENDIF
4404	ENDDO
4405	!
4406	! This is only to allow for the last check
4407	!
4408	inf = inflow_number(ing, inb)
4409	IF ( inflow_grid(ing, inb, inf) .EQ. ib .AND. inflow_basin(ing, inb, inf) .EQ. tokill) THEN
4410	doshift = .TRUE.
4411	ENDIF
4412	!
4413	IF ( .NOT. doshift ) THEN
4414	WRITE(numout,*) 'Strange we did not find the basin to kill in the downstream basin'
4415	STOP 'routing_killbas'
4416	ENDIF
4417	inflow_number(ing, inb) = inflow_number(ing, inb) - 1
4418	!
4419	ENDIF
4420	!
4421	! Now remove from the arrays the information of basin "tokill"
4422	!
4423	basin_id(ib, tokill:basin_count(ib)-1) = basin_id(ib, tokill+1:basin_count(ib))
4424	basin_flowdir(ib, tokill:basin_count(ib)-1) = basin_flowdir(ib, tokill+1:basin_count(ib))
4425	basin_area(ib, tokill:basin_count(ib)-1) = basin_area(ib, tokill+1:basin_count(ib))
4426	basin_area(ib, basin_count(ib):nwbas) = zero
4427	basin_topoind(ib, tokill:basin_count(ib)-1) = basin_topoind(ib, tokill+1:basin_count(ib))
4428	basin_topoind(ib, basin_count(ib):nwbas) = zero
4429	fetch_basin(ib, tokill:basin_count(ib)-1) = fetch_basin(ib, tokill+1:basin_count(ib))
4430	fetch_basin(ib, basin_count(ib):nwbas) = zero
4431	!
4432	! Before we remove the information from the outflow fields we have to correct the corresponding inflow fields
4433	! of the grids into which the flow goes
4434	!
4435	DO ibs = tokill+1,basin_count(ib)
4436	ing = outflow_grid(ib, ibs)
4437	inb = outflow_basin(ib, ibs)
4438	IF ( ing .GT. 0 ) THEN
4439	DO inf = 1, inflow_number(ing, inb)
4440	IF ( inflow_grid(ing,inb,inf) .EQ. ib .AND. inflow_basin(ing,inb,inf) .EQ. ibs) THEN
4441	inflow_basin(ing,inb,inf) = ibs - 1
4442	ENDIF
4443	ENDDO
4444	ENDIF
4445	ENDDO
4446	outflow_grid(ib, tokill:basin_count(ib)-1) = outflow_grid(ib, tokill+1:basin_count(ib))
4447	outflow_basin(ib, tokill:basin_count(ib)-1) = outflow_basin(ib, tokill+1:basin_count(ib))
4448	!
4449	! Basins which moved down also need to redirect their incoming flows.
4450	!
4451	DO ibs=tokill+1, basin_count(ib)
4452	DO inf = 1, inflow_number(ib, ibs)
4453	outflow_basin(inflow_grid(ib, ibs, inf), inflow_basin(ib, ibs, inf)) = ibs-1
4454	ENDDO
4455	ENDDO
4456	!
4457	! Shift the inflow basins
4458	!
4459	DO it = tokill+1,basin_count(ib)
4460	inflow_grid(ib, it-1, 1:inflow_number(ib,it)) = inflow_grid(ib, it, 1:inflow_number(ib,it))
4461	inflow_basin(ib, it-1, 1:inflow_number(ib,it)) = inflow_basin(ib, it, 1:inflow_number(ib,it))
4462	inflow_number(ib,it-1) = inflow_number(ib,it)
4463	ENDDO
4464	!
4465	basin_count(ib) = basin_count(ib) - 1
4466	!
4467	END SUBROUTINE routing_killbas
4468	!
4469	!---------------------------------------------------------------------------------
4470	!
4471	SUBROUTINE routing_names(numlar, basin_names)
4472	!
4473	IMPLICIT NONE
4474	!
4475	! Arguments
4476	!
4477	INTEGER(i_std), INTENT(in) :: numlar
4478	CHARACTER(LEN=*), INTENT(inout) :: basin_names(numlar)
4479	!
4480	! Local
4481	!
4482	INTEGER(i_std), PARAMETER :: listleng=349
4483	INTEGER(i_std) :: lenstr, i
4484	CHARACTER(LEN=60), DIMENSION(listleng) :: list_names
4485	CHARACTER(LEN=60) :: tmp_str
4486	!
4487	lenstr = LEN(basin_names(1))
4488	!
4489	list_names(1) = "Amazon"
4490	list_names(2) = "Nile"
4491	list_names(3) = "Zaire"
4492	list_names(4) = "Mississippi"
4493	list_names(5) = "Amur"
4494	list_names(6) = "Parana"
4495	list_names(7) = "Yenisei"
4496	list_names(8) = "Ob"
4497	list_names(9) = "Lena"
4498	list_names(10) = "Niger"
4499	list_names(11) = "Zambezi"
4500	list_names(12) = "Yangtze"
4501	list_names(13) = "Chang Jiang"
4502	list_names(14) = "Mackenzie"
4503	list_names(15) = "Ganges"
4504	list_names(16) = "Chari"
4505	list_names(17) = "Volga"
4506	list_names(18) = "St. Lawrence"
4507	list_names(19) = "Indus"
4508	list_names(20) = "Syr-Darya"
4509	list_names(21) = "Nelson"
4510	list_names(22) = "Orinoco"
4511	list_names(23) = "Murray"
4512	list_names(24) = "Great Artesian Basin"
4513	list_names(25) = "Shatt el Arab"
4514	list_names(26) = "Orange"
4515	list_names(27) = "Huang He"
4516	list_names(28) = "Yukon"
4517	list_names(29) = "Senegal"
4518	list_names(30) = "Chott Jerid"
4519	list_names(31) = "Jubba"
4520	list_names(32) = "Colorado (Ari)"
4521	list_names(33) = "Rio Grande (US)"
4522	list_names(34) = "Danube"
4523	list_names(35) = "Mekong"
4524	list_names(36) = "Tocantins"
4525	list_names(37) = "Wadi al Farigh"
4526	list_names(38) = "Tarim"
4527	list_names(39) = "Columbia"
4528	list_names(40) = "Noname (GHAASBasin49)"
4529	list_names(41) = "Kolyma"
4530	list_names(42) = "Sao Francisco"
4531	list_names(43) = "Amu-Darya"
4532	list_names(44) = "Noname (GHAASBasin51)"
4533	list_names(45) = "Dnepr"
4534	list_names(46) = "Noname (GHAASBasin61)"
4535	list_names(47) = "Don"
4536	list_names(48) = "Colorado (Arg)"
4537	list_names(49) = "Limpopo"
4538	list_names(50) = "GHAASBasin50"
4539	list_names(51) = "Zhujiang"
4540	list_names(52) = "Irrawaddy"
4541	list_names(53) = "Volta"
4542	list_names(54) = "GHAASBasin54"
4543	list_names(55) = "Farah"
4544	list_names(56) = "Khatanga"
4545	list_names(57) = "Dvina"
4546	list_names(58) = "Urugay"
4547	list_names(59) = "Qarqan"
4548	list_names(60) = "Noname (GHAASBasin75)"
4549	list_names(61) = "Parnaiba"
4550	list_names(62) = "Noname (GHAASBasin73)"
4551	list_names(63) = "Indigirka"
4552	list_names(64) = "Churchill (Hud)"
4553	list_names(65) = "Godavari"
4554	list_names(66) = "Pur - Taz"
4555	list_names(67) = "Pechora"
4556	list_names(68) = "Baker"
4557	list_names(69) = "Ural"
4558	list_names(70) = "Neva"
4559	list_names(71) = "Liao"
4560	list_names(72) = "Salween"
4561	list_names(73) = "GHAASBasin73"
4562	list_names(74) = "Jordan"
4563	list_names(75) = "Noname (GHAASBasin78)"
4564	list_names(76) = "Magdalena"
4565	list_names(77) = "Krishna"
4566	list_names(78) = "Salado"
4567	list_names(79) = "Fraser"
4568	list_names(80) = "Hai Ho"
4569	list_names(81) = "Huai"
4570	list_names(82) = "Yana"
4571	list_names(83) = "Noname (GHAASBasin95)"
4572	list_names(84) = "Noname (GHAASBasin105)"
4573	list_names(85) = "Kura"
4574	list_names(86) = "Olenek"
4575	list_names(87) = "Ogooue"
4576	list_names(88) = "Taymyr"
4577	list_names(89) = "Negro Arg"
4578	list_names(90) = "Chubut"
4579	list_names(91) = "GHAASBasin91"
4580	list_names(92) = "Noname (GHAASBasin122)"
4581	list_names(93) = "Noname (GHAASBasin120)"
4582	list_names(94) = "Sacramento"
4583	list_names(95) = "Fitzroy West"
4584	list_names(96) = "Grande de Santiago"
4585	list_names(97) = "Rufiji"
4586	list_names(98) = "Wisla"
4587	list_names(99) = "Noname (GHAASBasin47)"
4588	list_names(100) = "Noname (GHAASBasin127)"
4589	list_names(101) = "Hong"
4590	list_names(102) = "Noname (GHAASBasin97)"
4591	list_names(103) = "Swan-Avon"
4592	list_names(104) = "Rhine"
4593	list_names(105) = "Cuanza"
4594	list_names(106) = "GHAASBasin106"
4595	list_names(107) = "Noname (GHAASBasin142)"
4596	list_names(108) = "Roviuna"
4597	list_names(109) = "Essequibo"
4598	list_names(110) = "Elbe"
4599	list_names(111) = "Koksoak"
4600	list_names(112) = "Chao Phraya"
4601	list_names(113) = "Brahmani"
4602	list_names(114) = "Noname (GHAASBasin165)"
4603	list_names(115) = "Pyasina"
4604	list_names(116) = "Fitzroy East"
4605	list_names(117) = "Noname (GHAASBasin173)"
4606	list_names(118) = "Albany"
4607	list_names(119) = "Sanaga"
4608	list_names(120) = "GHAASBasin120"
4609	list_names(121) = "Noname (GHAASBasin178)"
4610	list_names(122) = "Noname (GHAASBasin148)"
4611	list_names(123) = "Brazos (Tex)"
4612	list_names(124) = "GHAASBasin124"
4613	list_names(125) = "Alabama"
4614	list_names(126) = "Noname (GHAASBasin174)"
4615	list_names(127) = "Noname (GHAASBasin179)"
4616	list_names(128) = "Balsas"
4617	list_names(129) = "Noname (GHAASBasin172)"
4618	list_names(130) = "Burdekin"
4619	list_names(131) = "Colorado (Texas)"
4620	list_names(132) = "Noname (GHAASBasin150)"
4621	list_names(133) = "Odra"
4622	list_names(134) = "Loire"
4623	list_names(135) = "Noname (GHAASBasin98)"
4624	list_names(136) = "Galana"
4625	list_names(137) = "Kuskowin"
4626	list_names(138) = "Moose"
4627	list_names(139) = "Narmada"
4628	list_names(140) = "GHAASBasin140"
4629	list_names(141) = "GHAASBasin141"
4630	list_names(142) = "Flinders"
4631	list_names(143) = "Kizil Irmak"
4632	list_names(144) = "GHAASBasin144"
4633	list_names(145) = "Save"
4634	list_names(146) = "Roper"
4635	list_names(147) = "Churchill (Atlantic)"
4636	list_names(148) = "GHAASBasin148"
4637	list_names(149) = "Victoria"
4638	list_names(150) = "Back"
4639	list_names(151) = "Bandama"
4640	list_names(152) = "Severn (Can)"
4641	list_names(153) = "Po"
4642	list_names(154) = "GHAASBasin154"
4643	list_names(155) = "GHAASBasin155"
4644	list_names(156) = "GHAASBasin156"
4645	list_names(157) = "Rhone"
4646	list_names(158) = "Tana (Ken)"
4647	list_names(159) = "La Grande"
4648	list_names(160) = "GHAASBasin160"
4649	list_names(161) = "Cunene"
4650	list_names(162) = "Douro"
4651	list_names(163) = "GHAASBasin163"
4652	list_names(164) = "Nemanus"
4653	list_names(165) = "GHAASBasin165"
4654	list_names(166) = "Anabar"
4655	list_names(167) = "Hayes"
4656	list_names(168) = "Mearim"
4657	list_names(169) = "GHAASBasin169"
4658	list_names(170) = "Panuco"
4659	list_names(171) = "GHAASBasin171"
4660	list_names(172) = "Doce"
4661	list_names(173) = "Gasgoyne"
4662	list_names(174) = "GHAASBasin174"
4663	list_names(175) = "GHAASBasin175"
4664	list_names(176) = "Ashburton"
4665	list_names(177) = "GHAASBasin177"
4666	list_names(178) = "Peel"
4667	list_names(179) = "Daugava"
4668	list_names(180) = "GHAASBasin180"
4669	list_names(181) = "Ebro"
4670	list_names(182) = "Comoe"
4671	list_names(183) = "Jacui"
4672	list_names(184) = "GHAASBasin184"
4673	list_names(185) = "Kapuas"
4674	list_names(186) = "GHAASBasin186"
4675	list_names(187) = "Penzhina"
4676	list_names(188) = "Cauweri"
4677	list_names(189) = "GHAASBasin189"
4678	list_names(190) = "Mamberamo"
4679	list_names(191) = "Sepik"
4680	list_names(192) = "GHAASBasin192"
4681	list_names(193) = "Sassandra"
4682	list_names(194) = "GHAASBasin194"
4683	list_names(195) = "GHAASBasin195"
4684	list_names(196) = "Nottaway"
4685	list_names(197) = "Barito"
4686	list_names(198) = "GHAASBasin198"
4687	list_names(199) = "Seine"
4688	list_names(200) = "Tejo"
4689	list_names(201) = "GHAASBasin201"
4690	list_names(202) = "Gambia"
4691	list_names(203) = "Susquehanna"
4692	list_names(204) = "Dnestr"
4693	list_names(205) = "Murchinson"
4694	list_names(206) = "Deseado"
4695	list_names(207) = "Mitchell"
4696	list_names(208) = "Mahakam"
4697	list_names(209) = "GHAASBasin209"
4698	list_names(210) = "Pangani"
4699	list_names(211) = "GHAASBasin211"
4700	list_names(212) = "GHAASBasin212"
4701	list_names(213) = "GHAASBasin213"
4702	list_names(214) = "GHAASBasin214"
4703	list_names(215) = "GHAASBasin215"
4704	list_names(216) = "Bug"
4705	list_names(217) = "GHAASBasin217"
4706	list_names(218) = "Usumacinta"
4707	list_names(219) = "Jequitinhonha"
4708	list_names(220) = "GHAASBasin220"
4709	list_names(221) = "Corantijn"
4710	list_names(222) = "Fuchun Jiang"
4711	list_names(223) = "Copper"
4712	list_names(224) = "Tapti"
4713	list_names(225) = "Menjiang"
4714	list_names(226) = "Karun"
4715	list_names(227) = "Mezen"
4716	list_names(228) = "Guadiana"
4717	list_names(229) = "Maroni"
4718	list_names(230) = "GHAASBasin230"
4719	list_names(231) = "Uda"
4720	list_names(232) = "GHAASBasin232"
4721	list_names(233) = "Kuban"
4722	list_names(234) = "Colville"
4723	list_names(235) = "Thaane"
4724	list_names(236) = "Alazeya"
4725	list_names(237) = "Paraiba do Sul"
4726	list_names(238) = "GHAASBasin238"
4727	list_names(239) = "Fortesque"
4728	list_names(240) = "GHAASBasin240"
4729	list_names(241) = "GHAASBasin241"
4730	list_names(242) = "Winisk"
4731	list_names(243) = "GHAASBasin243"
4732	list_names(244) = "GHAASBasin244"
4733	list_names(245) = "Ikopa"
4734	list_names(246) = "Gilbert"
4735	list_names(247) = "Kouilou"
4736	list_names(248) = "Fly"
4737	list_names(249) = "GHAASBasin249"
4738	list_names(250) = "GHAASBasin250"
4739	list_names(251) = "GHAASBasin251"
4740	list_names(252) = "Mangoky"
4741	list_names(253) = "Damodar"
4742	list_names(254) = "Onega"
4743	list_names(255) = "Moulouya"
4744	list_names(256) = "GHAASBasin256"
4745	list_names(257) = "Ord"
4746	list_names(258) = "GHAASBasin258"
4747	list_names(259) = "GHAASBasin259"
4748	list_names(260) = "GHAASBasin260"
4749	list_names(261) = "GHAASBasin261"
4750	list_names(262) = "Narva"
4751	list_names(263) = "GHAASBasin263"
4752	list_names(264) = "Seal"
4753	list_names(265) = "Cheliff"
4754	list_names(266) = "Garonne"
4755	list_names(267) = "Rupert"
4756	list_names(268) = "GHAASBasin268"
4757	list_names(269) = "Brahmani"
4758	list_names(270) = "Sakarya"
4759	list_names(271) = "Gourits"
4760	list_names(272) = "Sittang"
4761	list_names(273) = "Rajang"
4762	list_names(274) = "Evros"
4763	list_names(275) = "Appalachicola"
4764	list_names(276) = "Attawapiskat"
4765	list_names(277) = "Lurio"
4766	list_names(278) = "Daly"
4767	list_names(279) = "Penner"
4768	list_names(280) = "GHAASBasin280"
4769	list_names(281) = "GHAASBasin281"
4770	list_names(282) = "Guadalquivir"
4771	list_names(283) = "Nadym"
4772	list_names(284) = "GHAASBasin284"
4773	list_names(285) = "Saint John"
4774	list_names(286) = "GHAASBasin286"
4775	list_names(287) = "Cross"
4776	list_names(288) = "Omoloy"
4777	list_names(289) = "Oueme"
4778	list_names(290) = "GHAASBasin290"
4779	list_names(291) = "Gota"
4780	list_names(292) = "Nueces"
4781	list_names(293) = "Stikine"
4782	list_names(294) = "Yalu"
4783	list_names(295) = "Arnaud"
4784	list_names(296) = "GHAASBasin296"
4785	list_names(297) = "Jequitinhonha"
4786	list_names(298) = "Kamchatka"
4787	list_names(299) = "GHAASBasin299"
4788	list_names(300) = "Grijalva"
4789	list_names(301) = "GHAASBasin301"
4790	list_names(302) = "Kemijoki"
4791	list_names(303) = "Olifants"
4792	list_names(304) = "GHAASBasin304"
4793	list_names(305) = "Tsiribihina"
4794	list_names(306) = "Coppermine"
4795	list_names(307) = "GHAASBasin307"
4796	list_names(308) = "GHAASBasin308"
4797	list_names(309) = "Kovda"
4798	list_names(310) = "Trinity"
4799	list_names(311) = "Glama"
4800	list_names(312) = "GHAASBasin312"
4801	list_names(313) = "Luan"
4802	list_names(314) = "Leichhardt"
4803	list_names(315) = "GHAASBasin315"
4804	list_names(316) = "Gurupi"
4805	list_names(317) = "GR Baleine"
4806	list_names(318) = "Aux Feuilles"
4807	list_names(319) = "GHAASBasin319"
4808	list_names(320) = "Weser"
4809	list_names(321) = "GHAASBasin321"
4810	list_names(322) = "GHAASBasin322"
4811	list_names(323) = "Yesil"
4812	list_names(324) = "Incomati"
4813	list_names(325) = "GHAASBasin325"
4814	list_names(326) = "GHAASBasin326"
4815	list_names(327) = "Pungoe"
4816	list_names(328) = "GHAASBasin328"
4817	list_names(329) = "Meuse"
4818	list_names(330) = "Eastmain"
4819	list_names(331) = "Araguari"
4820	list_names(332) = "Hudson"
4821	list_names(333) = "GHAASBasin333"
4822	list_names(334) = "GHAASBasin334"
4823	list_names(335) = "GHAASBasin335"
4824	list_names(336) = "GHAASBasin336"
4825	list_names(337) = "Kobuk"
4826	list_names(338) = "Altamaha"
4827	list_names(339) = "GHAASBasin339"
4828	list_names(340) = "Mand"
4829	list_names(341) = "Santee"
4830	list_names(342) = "GHAASBasin342"
4831	list_names(343) = "GHAASBasin343"
4832	list_names(344) = "GHAASBasin344"
4833	list_names(345) = "Hari"
4834	list_names(346) = "GHAASBasin346"
4835	list_names(347) = "Wami"
4836	list_names(348) = "GHAASBasin348"
4837	list_names(349) = "GHAASBasin349"
4838	!
4839	basin_names(:) = ' '
4840	!
4841	DO i=1,numlar
4842	tmp_str = list_names(i)
4843	basin_names(i) = tmp_str(1:MIN(lenstr,LEN_TRIM(tmp_str)))
4844	ENDDO
4845	!
4846	END SUBROUTINE routing_names
4847	!
4848	!---------------------------------------------------------------------------------
4849	!
4850	SUBROUTINE routing_irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, &
4851	& init_irrig, irrigated, init_flood, floodplains, hist_id, hist2_id)
4852	!
4853	! This program will interpoalte the 0.5 x 05 deg based data set to the resolution
4854	! of the model.
4855	!
4856	IMPLICIT NONE
4857	!
4858	! INPUT
4859	!
4860	INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated
4861	INTEGER(i_std), INTENT(in) :: index(nbpt) ! Index on the global map.
4862	REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !)
4863	INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
4864	REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in m of each grid-box in X and Y
4865	REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box
4866	LOGICAL, INTENT(in) :: init_irrig ! Do we need to compute irrigation ?
4867	REAL(r_std), INTENT(out) :: irrigated(:) ! Irrigated surface in each grid-box
4868	LOGICAL, INTENT(in) :: init_flood ! Do we need to compute floodplains
4869	REAL(r_std), INTENT(out) :: floodplains(:) ! Surface which can be inondated in each grid-box
4870	INTEGER(i_std), INTENT(in) :: hist_id ! Access to history file
4871	INTEGER(i_std), INTENT(in) :: hist2_id ! Access to history file 2
4872	!
4873	! LOCAL
4874	!
4875	CHARACTER(LEN=80) :: filename
4876	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, ilf, lastjp, nbexp
4877	REAL(r_std) :: lev(1), date, dt, coslat
4878	REAL(r_std) :: lon_up, lon_low, lat_up, lat_low
4879	INTEGER(i_std) :: itau(1)
4880	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, lat_ful, lon_ful
4881	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrigated_frac, flood_frac
4882	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: loup_rel, lolow_rel, laup_rel, lalow_rel
4883	REAL(r_std) :: area_irrig, area_flood, ax, ay, sgn, surp
4884	REAL(r_std) :: lonrel, louprel, lolowrel
4885	REAL(r_std) :: irrigmap(nbpt)
4886	!
4887	!
4888	!Config Key = IRRIGATION_FILE
4889	!Config Desc = Name of file which contains the map of irrigated areas
4890	!Config Def = irrigated.nc
4891	!Config If = IRRIGATE
4892	!Config Help = The name of the file to be opened to read the field
4893	!Config with the area in m^2 of the area irrigated within each
4894	!Config 0.5 0.5 deg grid box. The map currently used is the one
4895	!Config developed by the Center for Environmental Systems Research
4896	!Config in Kassel (1995).
4897	!
4898	filename = 'irrigated.nc'
4899	CALL getin_p('IRRIGATION_FILE',filename)
4900	!
4901	IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid)
4902	CALL bcast(iml)
4903	CALL bcast(jml)
4904	CALL bcast(lml)
4905	CALL bcast(tml)
4906	!
4907	!
4908	ALLOCATE (lat_rel(iml,jml))
4909	ALLOCATE (lon_rel(iml,jml))
4910	ALLOCATE (laup_rel(iml,jml))
4911	ALLOCATE (loup_rel(iml,jml))
4912	ALLOCATE (lalow_rel(iml,jml))
4913	ALLOCATE (lolow_rel(iml,jml))
4914	ALLOCATE (lat_ful(iml+2,jml+2))
4915	ALLOCATE (lon_ful(iml+2,jml+2))
4916	ALLOCATE (irrigated_frac(iml,jml))
4917	ALLOCATE (flood_frac(iml,jml))
4918	!
4919	IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid)
4920	CALL bcast(lon_rel)
4921	CALL bcast(lat_rel)
4922	CALL bcast(lev)
4923	CALL bcast(itau)
4924	CALL bcast(date)
4925	CALL bcast(dt)
4926
4927	!
4928	IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 1, 1, irrigated_frac)
4929	CALL bcast(irrigated_frac)
4930	IF (is_root_prc) CALL flinget(fid, 'inundated', iml, jml, lml, tml, 1, 1, flood_frac)
4931	CALL bcast(flood_frac)
4932	!
4933	IF (is_root_prc) CALL flinclo(fid)
4934	!
4935	! Set to zero all fraction which are less than 0.5%
4936	!
4937	DO ip=1,iml
4938	DO jp=1,jml
4939	!
4940	IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN
4941	irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100.
4942	IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero
4943	ENDIF
4944	!
4945	IF ( flood_frac(ip,jp) .LT. undef_sechiba-un) THEN
4946	flood_frac(ip,jp) = flood_frac(ip,jp)/100
4947	IF ( flood_frac(ip,jp) < 0.005 ) flood_frac(ip,jp) = zero
4948	ENDIF
4949	!
4950	ENDDO
4951	ENDDO
4952	!
4953	WRITE(numout,*) 'lon_rel : ', MAXVAL(lon_rel), MINVAL(lon_rel)
4954	WRITE(numout,*) 'lat_rel : ', MAXVAL(lat_rel), MINVAL(lat_rel)
4955	WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), &
4956	& MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba)
4957	WRITE(numout,*) 'flood_frac : ', MINVAL(flood_frac, MASK=flood_frac .GT. 0), &
4958	& MAXVAL(flood_frac, MASK=flood_frac .LT. undef_sechiba)
4959	!
4960	nbexp = 0
4961	!
4962	! Duplicate the border assuming we have a global grid going from west to east
4963	!
4964	lon_ful(2:iml+1,2:jml+1) = lon_rel(1:iml,1:jml)
4965	lat_ful(2:iml+1,2:jml+1) = lat_rel(1:iml,1:jml)
4966	!
4967	IF ( lon_rel(iml,1) .LT. lon_ful(2,2)) THEN
4968	lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)
4969	lat_ful(1,2:jml+1) = lat_rel(iml,1:jml)
4970	ELSE
4971	lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)-360
4972	lat_ful(1,2:jml+1) = lat_rel(iml,1:jml)
4973	ENDIF
4974
4975	IF ( lon_rel(1,1) .GT. lon_ful(iml+1,2)) THEN
4976	lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)
4977	lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml)
4978	ELSE
4979	lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)+360
4980	lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml)
4981	ENDIF
4982	!
4983	sgn = lat_rel(1,1)/ABS(lat_rel(1,1))
4984	lat_ful(2:iml+1,1) = sgn*180 - lat_rel(1:iml,1)
4985	sgn = lat_rel(1,jml)/ABS(lat_rel(1,jml))
4986	lat_ful(2:iml+1,jml+2) = sgn*180 - lat_rel(1:iml,jml)
4987	lat_ful(1,1) = lat_ful(iml+1,1)
4988	lat_ful(iml+2,1) = lat_ful(2,1)
4989	lat_ful(1,jml+2) = lat_ful(iml+1,jml+2)
4990	lat_ful(iml+2,jml+2) = lat_ful(2,jml+2)
4991	!
4992	! Add the longitude lines to the top and bottom
4993	!
4994	lon_ful(:,1) = lon_ful(:,2)
4995	lon_ful(:,jml+2) = lon_ful(:,jml+1)
4996	!
4997	! Get the upper and lower limits of each grid box
4998	!
4999	DO ip=1,iml
5000	DO jp=1,jml
5001	loup_rel(ip,jp) =MAX(0.5(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
5002	lolow_rel(ip,jp) =MIN(0.5(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
5003	laup_rel(ip,jp) =MAX(0.5(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
5004	lalow_rel(ip,jp) =MIN(0.5(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
5005	ENDDO
5006	ENDDO
5007	!
5008	! Now we take each grid point and find out which values from the forcing we need to average
5009	!
5010	DO ib =1, nbpt
5011	!
5012	! We find the 4 limits of the grid-box. As we transform the resolution of the model
5013	! into longitudes and latitudes we do not have the problem of periodicity.
5014	! coslat is a help variable here !
5015	!
5016	coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )pi/180. R_Earth
5017	!
5018	lon_up = lalo(ib,2)+ resolution(ib,1)/(2.0*coslat)
5019	lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat)
5020	!
5021	coslat = pi/180. * R_Earth
5022	!
5023	lat_up =lalo(ib,1)+resolution(ib,2)/(2.0*coslat)
5024	lat_low =lalo(ib,1)-resolution(ib,2)/(2.0*coslat)
5025	!
5026	!
5027	! Find the grid boxes from the data that go into the model's boxes
5028	! We still work as if we had a regular grid ! Well it needs to be localy regular so
5029	! so that the longitude at the latitude of the last found point is close to the one of the next point.
5030	!
5031	lastjp = 1
5032	area_irrig = zero
5033	area_flood = zero
5034	DO ip=1,iml
5035	!
5036	! Either the center of the data grid point is in the interval of the model grid or
5037	! the East and West limits of the data grid point are on either sides of the border of
5038	! the data grid.
5039	!
5040	!
5041	! We find the 4 limits of the grid-box. As we transform the resolution of the model
5042	! into longitudes and latitudes we do not have the problem of periodicity.
5043	! coslat is a help variable here !
5044	!
5045	!
5046	! To do that correctly we have to check if the grid box sits on the date-line.
5047	!
5048	IF ( lon_low < -180.0 ) THEN
5049	lonrel = MOD( lon_rel(ip,lastjp) - 360.0, 360.0)
5050	lolowrel = MOD( lolow_rel(ip,lastjp) - 360.0, 360.0)
5051	louprel = MOD( loup_rel(ip,lastjp) - 360.0, 360.0)
5052	!
5053	ELSE IF ( lon_up > 180.0 ) THEN
5054	lonrel = MOD( 360. - lon_rel(ip,lastjp), 360.0)
5055	lolowrel = MOD( 360. - lolow_rel(ip,lastjp), 360.0)
5056	louprel = MOD( 360. - loup_rel(ip,lastjp), 360.0)
5057	ELSE
5058	lonrel = lon_rel(ip,lastjp)
5059	lolowrel = lolow_rel(ip,lastjp)
5060	louprel = loup_rel(ip,lastjp)
5061	ENDIF
5062	!
5063	!
5064	!
5065	IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
5066	& lolowrel < lon_low .AND. louprel > lon_low .OR. &
5067	& lolowrel < lon_up .AND. louprel > lon_up ) THEN
5068	!
5069	DO jp = 1, jml
5070	!
5071	! Now that we have the longitude let us find the latitude
5072	!
5073	IF ( lat_rel(ip,jp) > lat_low .AND. lat_rel(ip,jp) < lat_up .OR. &
5074	& lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.&
5075	& lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN
5076	!
5077	lastjp = jp
5078	!
5079	! Mising values in the file are assumed to be 1e20
5080	!
5081	IF ( lon_low < -180.0 ) THEN
5082	lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
5083	louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
5084	!
5085	ELSE IF ( lon_up > 180.0 ) THEN
5086	lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0)
5087	louprel = MOD( 360. - loup_rel(ip,jp), 360.0)
5088	ELSE
5089	lolowrel = lolow_rel(ip,jp)
5090	louprel = loup_rel(ip,jp)
5091	ENDIF
5092	!
5093	! Get the area of the fine grid in the model grid
5094	!
5095	coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), 0.001 )
5096	ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))pi/180. R_Earth * coslat
5097	ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))pi/180. R_Earth
5098	!
5099	IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN
5100	area_irrig = area_irrig + axayirrigated_frac(ip,jp)
5101	ENDIF
5102	!
5103	IF (flood_frac(ip,jp) .LT. undef_sechiba-un) THEN
5104	area_flood = area_flood + axayflood_frac(ip,jp)
5105	ENDIF
5106	!
5107	ENDIF
5108	!
5109	ENDDO
5110	!
5111	ENDIF
5112	!
5113	ENDDO
5114	!
5115	! Put the toal irrigated and flooded areas in the output variables
5116	!
5117	IF ( init_irrig ) THEN
5118	irrigated(ib) = MIN(area_irrig, resolution(ib,1)resolution(ib,2)contfrac(ib))
5119	IF ( irrigated(ib) < 0 ) THEN
5120	WRITE(numout,*) 'We have a problem here : ', irrigated(ib)
5121	WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
5122	WRITE(numout,*) area_irrig
5123	STOP
5124	ENDIF
5125	! Compute a diagnostic of the map.
5126	IF(contfrac(ib).GT.zero) THEN
5127	irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)resolution(ib,2)contfrac(ib) )
5128	ELSE
5129	irrigmap (ib) = zero
5130	ENDIF
5131	!
5132	ENDIF
5133	!
5134	IF ( init_flood ) THEN
5135	floodplains(ib) = MIN(area_flood, resolution(ib,1)resolution(ib,2)contfrac(ib))
5136	IF ( floodplains(ib) < 0 ) THEN
5137	WRITE(numout,*) 'We have a problem here : ', floodplains(ib)
5138	WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2)
5139	WRITE(numout,*) area_flood
5140	STOP
5141	ENDIF
5142	ENDIF
5143	!
5144	!
5145	ENDDO
5146	!
5147	IF ( init_irrig .AND. init_flood ) THEN
5148	DO ib = 1, nbpt
5149	IF ( floodplains(ib)+irrigated(ib) > resolution(ib,1)resolution(ib,2)contfrac(ib) ) THEN
5150	surp = (resolution(ib,1)resolution(ib,2)contfrac(ib))/(floodplains(ib)+irrigated(ib))
5151	floodplains(ib) = floodplains(ib) * surp
5152	irrigated(ib) = irrigated(ib) * surp
5153	ENDIF
5154	ENDDO
5155	ENDIF
5156	!
5157	IF ( init_irrig ) THEN
5158	WRITE(numout,*) 'RESULT irrigated : ', MINVAL(irrigated), MAXVAL(irrigated)
5159	IF ( .NOT. almaoutput ) THEN
5160	CALL histwrite(hist_id, 'irrigmap', 1, irrigmap, nbpt, index)
5161	ELSE
5162	CALL histwrite(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index)
5163	ENDIF
5164	IF ( hist2_id > 0 ) THEN
5165	IF ( .NOT. almaoutput ) THEN
5166	CALL histwrite(hist2_id, 'irrigmap', 1, irrigmap, nbpt, index)
5167	ELSE
5168	CALL histwrite(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index)
5169	ENDIF
5170	ENDIF
5171	ENDIF
5172	!
5173	IF ( init_flood ) THEN
5174	WRITE(numout,*) 'RESULT floodplains : ', MINVAL(floodplains), MAXVAL(floodplains)
5175	ENDIF
5176	!
5177	END SUBROUTINE routing_irrigmap
5178	!
5179	!---------------------------------------------------------------------------------
5180	!
5181	SUBROUTINE routing_waterbal(nbpt, firstcall, runoff, drainage, returnflow, &
5182	& irrigation, riverflow, coastalflow)
5183	!
5184	IMPLICIT NONE
5185	!
5186	! This subroutine should allow us to check the water balance in the routing module.
5187	!
5188	INTEGER(i_std), INTENT(in) :: nbpt ! Domain size
5189	LOGICAL, INTENT(in) :: firstcall ! controls behaviour
5190	REAL(r_std), INTENT(in) :: runoff(nbpt) ! grid-point runoff
5191	REAL(r_std), INTENT(in) :: drainage(nbpt) ! grid-point drainage
5192	REAL(r_std), INTENT(in) :: returnflow(nbpt) ! The water flow which returns to the grid box (kg/m^2/dt)
5193	REAL(r_std), INTENT(in) :: irrigation(nbpt) ! Irrigation flux (kg/m^2 per dt)
5194	REAL(r_std), INTENT(in) :: riverflow(nbpt) ! Outflow of the major rivers (kg/dt)
5195	REAL(r_std), INTENT(in) :: coastalflow(nbpt) ! Outflow on coastal points by small basins (kg/dt)
5196	!
5197	! We sum-up all the water we have in the warious reservoirs
5198	!
5199	REAL(r_std), SAVE :: totw_stream, totw_fast, totw_slow, totw_lake
5200	REAL(r_std), SAVE :: totw_in, totw_out, totw_return, totw_irrig, totw_river, totw_coastal
5201	REAL(r_std) :: totarea, area
5202	!
5203	! Just to make sure we do not get too large numbers !
5204	!
5205	REAL(r_std), PARAMETER :: scaling = 1.0E+6
5206	REAL(r_std), PARAMETER :: allowed_err = 50.
5207	!
5208	INTEGER(i_std) :: ig
5209	!
5210	IF ( firstcall ) THEN
5211	!
5212	totw_stream = zero
5213	totw_fast = zero
5214	totw_slow = zero
5215	totw_lake = zero
5216	totw_in = zero
5217	!
5218	DO ig=1,nbpt
5219	!
5220	totarea = SUM(routing_area(ig,:))
5221	!
5222	totw_stream = totw_stream + SUM(stream_reservoir(ig,:)/scaling)
5223	totw_fast = totw_fast + SUM(fast_reservoir(ig,:)/scaling)
5224	totw_slow = totw_slow + SUM(slow_reservoir(ig,:)/scaling)
5225	totw_lake = totw_lake + lake_reservoir(ig)/scaling
5226	!
5227	totw_in = totw_in + (runoff(ig)totarea + drainage(ig)totarea)/scaling
5228	!
5229	ENDDO
5230	!
5231	ELSE
5232	!
5233	totw_out = zero
5234	totw_return = zero
5235	totw_irrig = zero
5236	totw_river = zero
5237	totw_coastal = zero
5238	area = zero
5239	!
5240	DO ig=1,nbpt
5241	!
5242	totarea = SUM(routing_area(ig,:))
5243	!
5244	totw_stream = totw_stream - SUM(stream_reservoir(ig,:)/scaling)
5245	totw_fast = totw_fast - SUM(fast_reservoir(ig,:)/scaling)
5246	totw_slow = totw_slow - SUM(slow_reservoir(ig,:)/scaling)
5247	totw_lake = totw_lake - lake_reservoir(ig)/scaling
5248	!
5249	totw_return = totw_return + returnflow(ig)*totarea/scaling
5250	totw_irrig = totw_irrig + irrigation(ig)*totarea/scaling
5251	totw_river = totw_river + riverflow(ig)/scaling
5252	totw_coastal = totw_coastal + coastalflow(ig)/scaling
5253	!
5254	area = area + totarea
5255	!
5256	ENDDO
5257	totw_out = totw_return + totw_irrig + totw_river + totw_coastal
5258	!
5259	! Now we have all the information to balance our water
5260	!
5261	IF ( ABS((totw_stream + totw_fast + totw_slow + totw_lake) - (totw_out - totw_in)) > allowed_err ) THEN
5262	WRITE(numout,*) 'WARNING : Water not conserved in routing. Limit at ', allowed_err, ' 10^6 kg'
5263	WRITE(numout,*) '--Water-- change : stream fast ', totw_stream, totw_fast
5264	WRITE(numout,*) '--Water-- change : slow, lake ', totw_slow, totw_lake
5265	WRITE(numout,*) '--Water>>> change in the routing res. : ', totw_stream + totw_fast + totw_slow + totw_lake
5266	WRITE(numout,*) '--Water input : ', totw_in
5267	WRITE(numout,*) '--Water output : ', totw_out
5268	WRITE(numout,*) '--Water output : return, irrig ', totw_return, totw_irrig
5269	WRITE(numout,*) '--Water output : river, coastal ',totw_river, totw_coastal
5270	WRITE(numout,*) '--Water>>> change by fluxes : ', totw_out - totw_in, ' Diff [mm/dt]: ', &
5271	& ((totw_stream + totw_fast + totw_slow + totw_lake) - (totw_out - totw_in))/area
5272	ENDIF
5273	!
5274	ENDIF
5275	!
5276	END SUBROUTINE routing_waterbal
5277	!
5278	!
5279	END MODULE routing

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

Download in other formats: