Context Navigation

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

Last change on this file since 282 was 277, checked in by didier.solyga, 13 years ago
Update the externalized version with the last commit of the trunk (revision 275)
File size: 209.4 KB

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

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

Download in other formats: