Context Navigation

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

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

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

Download in other formats: