Context Navigation

source: perso/abdelouhab.djerrah/ORCHIDEE/src_sechiba/routing.f90 @ 938

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

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

Download in other formats: