Context Navigation

source: tags/ORCHIDEE/src_sechiba/routing.f90 @ 6

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

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

Download in other formats: