Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/routing_tools.f90 @ 5143

Last change on this file since 5143 was 3845, checked in by trung.nguyen, 8 years ago
Corrected a tiny bug in routing_reg_bestsubbasin.
File size: 35.0 KB

Line
1	! =================================================================================================================================
2	! MODULE : routing_tools
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2016)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF
10	!!
11	!!\n DESCRIPTION: None
12	!!
13	!! RECENT CHANGE(S): None
14	!!
15	!! REFERENCE(S) :
16	!!
17	!! SVN :
18	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-ROUTING/ORCHIDEE/src_sechiba/routing.f90 $
19	!! $Date: 2016-05-19 15:16:26 +0200 (Thu, 19 May 2016) $
20	!! $Revision: 3449 $
21	!! \n
22	!_ ================================================================================================================================
23	!
24	MODULE routing_tools
25
26	USE ioipsl
27	USE xios_orchidee
28	USE ioipsl_para
29	USE constantes
30	USE grid
31	USE mod_orchidee_para
32
33	IMPLICIT NONE
34	PUBLIC :: routing_nextgrid_g, routing_anglediff_g
35	! We are working on the rectengular lat/lon grid of routing.nc
36	! The basic information of these grids are saved here.
37	!
38	INTEGER(i_std), PARAMETER :: nbne=8
39	! The hypothesis of the flow directions used in routing.nc :
40	REAL(r_std), PARAMETER, DIMENSION(nbne) :: rose=(/0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0/)
41	!
42	CONTAINS
43	!! ================================================================================================================================
44	!! SUBROUTINE : routing_nextgrid_g
45	!!
46	!>\BRIEF This function finds the grid box into which the flow direction points and we should
47	!! be able to find the next sub-basin of the river.
48	!! We will assume that this function is called on the root processor and with the global index of the
49	!! grid point.
50	!!
51	!! DESCRIPTION (definitions, functional, design, flags) : None
52	!!
53	!! RECENT CHANGE(S): None
54	!!
55	!! MAIN OUTPUT VARIABLE(S):
56	!!
57	!! REFERENCES : None
58	!!
59	!! FLOWCHART : None
60	!! \n
61	!_ ================================================================================================================================
62	!
63	INTEGER(i_std) FUNCTION routing_nextgrid_g(pnt, trip, inc)
64	!
65	! Input
66	!
67	INTEGER(i_std) :: pnt, trip
68	INTEGER(i_std), OPTIONAL :: inc
69	!
70	! Local
71	!
72	INTEGER(i_std) :: i
73	REAL(r_std), DIMENSION(NbNeighb) :: diffangle
74	INTEGER(i_std), DIMENSION(1) :: imin
75	INTEGER(i_std) :: trip_loc, nbind
76	!
77	! Compute the angle between the direction of the routing.nc flow and the headings perpendicular of the segments
78	! of the polygon.
79	!
80	IF (.NOT. is_root_prc) THEN
81	WRITE(,) "The function routing_nextgrid_g is written to work on the root proc and the full grid"
82	STOP "routing_nextgrid_g"
83	ENDIF
84	!
85	! Make sure that the flow direction coming from routing.f90 is between 1 and 8. This allows to call the
86	! function without having to make sure trip is between 1 and 8.
87	!
88	trip_loc = MODULO(trip-1, nbne)+1
89	!
90	DO i=1,NbNeighb
91	diffangle(i) = ABS(haversine_diffangle(headings_g(pnt,i), rose(trip_loc)))
92	ENDDO
93	!
94	imin = MINLOC(diffangle)
95	!
96	IF ( PRESENT(inc) ) THEN
97	! Ensure that the increment the users asks is within the NbNeighb the polygone has.
98	nbind = MODULO(imin(1)+inc-1, NbNeighb)+1
99	ELSE
100	nbind = imin(1)
101	ENDIF
102	!
103	routing_nextgrid_g = neighbours_g(pnt,nbind)
104	!
105	END FUNCTION routing_nextgrid_g
106	!
107	!! ================================================================================================================================
108	!! SUBROUTINE : routing_nextgrid
109	!!
110	!>\BRIEF This function finds the grid box into which the flow direction points and we should
111	!! be able to find the next sub-basin of the river.
112	!! We will assume that this function is called on the root processor and with the global index of the
113	!! grid point. This verison of the function can be called in a parallel region.
114	!!
115	!! DESCRIPTION (definitions, functional, design, flags) : None
116	!!
117	!! RECENT CHANGE(S): None
118	!!
119	!! MAIN OUTPUT VARIABLE(S):
120	!!
121	!! REFERENCES : None
122	!!
123	!! FLOWCHART : None
124	!! \n
125	!_ ================================================================================================================================
126	!
127	INTEGER(i_std) FUNCTION routing_nextgrid(pnt, trip, inc)
128	!
129	! Input
130	!
131	INTEGER(i_std) :: pnt, trip
132	INTEGER(i_std), OPTIONAL :: inc
133	!
134	! Local
135	!
136	INTEGER(i_std) :: i
137	REAL(r_std), DIMENSION(NbNeighb) :: diffangle
138	INTEGER(i_std), DIMENSION(1) :: imin
139	INTEGER(i_std) :: trip_loc, nbind
140	!
141	! Make sure that the flow direction coming from routing.f90 is between 1 and 8. This allows to call the
142	! function without having to make sure trip is between 1 and 8.
143	!
144	trip_loc = MODULO(trip-1, nbne)+1
145	!
146	! Compute the angle between the direction of the routing.nc flow and the headings perpendicular of the segments
147	! of the polygon.
148	!
149	DO i=1,NbNeighb
150	diffangle(i) = ABS(haversine_diffangle(headings(pnt,i), rose(trip_loc)))
151	ENDDO
152	!
153	imin = MINLOC(diffangle)
154	!
155	IF ( PRESENT(inc) ) THEN
156	! Ensure that the increment the users asks is within the NbNeighb the polygone has.
157	nbind = MODULO(imin(1)+inc-1, NbNeighb)+1
158	ELSE
159	nbind = imin(1)
160	ENDIF
161	!
162	routing_nextgrid = neighbours(pnt,nbind)
163	!
164	END FUNCTION routing_nextgrid
165	!
166	!
167	!! ================================================================================================================================
168	!! SUBROUTINE : routing_anglediff_g
169	!!
170	!>\BRIEF This routine computes the angle difference between the flow direction on the routing.nc grid (ytip) and
171	!! the heading of the line which links the 2 points (pnt and pntout) of the model grid. These 2 grid
172	!! points have to be neighbours.
173	!!
174	!! DESCRIPTION (definitions, functional, design, flags) : None
175	!!
176	!! RECENT CHANGE(S): None
177	!!
178	!! MAIN OUTPUT VARIABLE(S):
179	!!
180	!! REFERENCES : None
181	!!
182	!! FLOWCHART : None
183	!! \n
184	!_ ================================================================================================================================
185	!
186	REAL(r_std) FUNCTION routing_anglediff_g(pnt, pntout, trip)
187	!
188	! Input
189	!
190	INTEGER(i_std) :: pnt, pntout
191	INTEGER(i_std) :: trip
192	!
193	! Local
194	!
195	INTEGER(i_std) :: i
196	INTEGER(i_std) :: imin
197	INTEGER(i_std) :: trip_loc
198	!
199	! This version of the routine can only be called on the root proc.
200	!
201	IF (.NOT. is_root_prc) THEN
202	WRITE(,) "The function routing_flowangle_g is written to work on the root proc and the full grid"
203	STOP "routing_flowangle_g"
204	ENDIF
205	trip_loc = MOD(trip+nbne-1, nbne)+1
206	!
207	! Find the index of the neighbour we are considering so we get the heading out of pnt
208	!
209	imin=-1
210	DO i=1,NbNeighb
211	IF ( neighbours_g(pnt,i) .EQ. pntout ) THEN
212	imin = i
213	ENDIF
214	ENDDO
215	IF ( i > 0 ) THEN
216	routing_anglediff_g = ABS(haversine_diffangle(headings_g(pnt,imin), rose(trip_loc)))
217	ELSE
218	WRITE(,) "Input arguments : pnt, pntout, trip = ", pnt, pntout, trip
219	CALL ipslerr_p(3,'routing_anglediff_g', &
220	"Cannot find the neighbour requested.", &
221	& "None of the neighbours of the specified point (pnt) match pntout.", "")
222	ENDIF
223	!
224	END FUNCTION routing_anglediff_g
225	!
226	!
227	!! ================================================================================================================================
228	!! SUBROUTINE : routing_anglediff
229	!!
230	!>\BRIEF This routine computes the angle difference between the flow direction on the routing.nc grid (ytip) and
231	!! the heading of the line which links the 2 points (pnt and pntout) of the model grid. These 2 grid
232	!! points have to be neighbours.
233	!!
234	!! DESCRIPTION (definitions, functional, design, flags) : None
235	!!
236	!! RECENT CHANGE(S): None
237	!!
238	!! MAIN OUTPUT VARIABLE(S):
239	!!
240	!! REFERENCES : None
241	!!
242	!! FLOWCHART : None
243	!! \n
244	!_ ================================================================================================================================
245	!
246	REAL(r_std) FUNCTION routing_anglediff(pnt, pntout, trip)
247	!
248	! Input
249	!
250	INTEGER(i_std) :: pnt, pntout
251	INTEGER(i_std) :: trip
252	!
253	! Local
254	!
255	INTEGER(i_std) :: i
256	INTEGER(i_std) :: imin
257	INTEGER(i_std) :: trip_loc
258	!
259	! This version of the routine can only be called on the root proc.
260	!
261	trip_loc = MOD(trip+nbne-1, nbne)+1
262	!
263	! Find the index of the neighbour we are considering so we get the heading out of pnt
264	!
265	imin=-1
266	DO i=1,NbNeighb
267	IF ( neighbours(pnt,i) .EQ. pntout ) THEN
268	imin = i
269	ENDIF
270	ENDDO
271	IF ( i > 0 ) THEN
272	routing_anglediff = ABS(haversine_diffangle(headings(pnt,imin), rose(trip_loc)))
273	ELSE
274	WRITE(,) "Input arguments : pnt, pntout, trip = ", pnt, pntout, trip
275	CALL ipslerr_p(3,'routing_anglediff', &
276	"Cannot find the neighbour requested.", &
277	& "None of the neighbours of the specified point (pnt) match pntout.", "")
278	ENDIF
279	!
280	END FUNCTION routing_anglediff
281	!
282	!
283	!!
284	!! ================================================================================================================================
285	!! SUBROUTINE : routing_bestsubbasin()
286	!!
287	!>\BRIEF This routine help routing_linkup to find the most suitable sub-basin in the grid box
288	!! in which the source sub-basin flows.
289	!!
290	!! DESCRIPTION (definitions, functional, design, flags) : None
291	!!
292	!! RECENT CHANGE(S): None
293	!!
294	!! MAIN OUTPUT VARIABLE(S):
295	!!
296	!! REFERENCES : None
297	!!
298	!! FLOWCHART : None
299	!! \n
300	!_ ================================================================================================================================
301	!
302	SUBROUTINE routing_bestsubbasin(src_grid, src_id, src_hierarchy, src_flowdir, invented_basins, &
303	& ptmax, basmax, pt, nbsubbas, ids, hierarchy, flowdir, nbcoastals, coastals, &
304	& outbas, quality)
305	!
306	! Input
307	!
308	INTEGER(i_std), INTENT(in) :: src_grid, src_id, src_flowdir
309	REAL(r_std), INTENT(in) :: src_hierarchy
310	REAL(r_std), INTENT(in) :: invented_basins
311	!
312	INTEGER(i_std), INTENT(in) :: ptmax, basmax, pt
313	INTEGER(i_std), INTENT(in) :: nbsubbas(ptmax)
314	INTEGER(i_std), INTENT(in) :: ids(ptmax,basmax), flowdir(ptmax,basmax)
315	REAL(r_std), INTENT(in) :: hierarchy(ptmax,basmax)
316	INTEGER(i_std), INTENT(in) :: nbcoastals(ptmax)
317	INTEGER(i_std), INTENT(in) :: coastals(ptmax,basmax)
318	!
319	! Output
320	!
321	INTEGER(i_std), INTENT(out) :: outbas
322	REAL(r_std), INTENT(out) :: quality
323	!
324	! Local
325	!
326	INTEGER(i_std) :: i, sbl, nbfbas, nbopt
327	INTEGER(i_std), DIMENSION(1) :: ff
328	INTEGER(i_std), DIMENSION(basmax) :: fbas, tbas, sboptions, angle, subbas
329	INTEGER(i_std), DIMENSION(basmax) :: fbas_flowdir, tbas_flowdir
330	REAL(r_std), DIMENSION(basmax) :: fbas_hierarchy, tbas_hierarchy
331	!
332	! Set some default values
333	!
334	nbfbas = 0
335	nbopt = 0
336	sboptions(:) = undef_int
337	outbas = undef_int
338	!
339	! 1.0 Find possible basins in the list of the target grid
340	!
341	DO sbl=1,nbsubbas(pt)
342	!
343	! Either it is a standard basin or one aggregated from ocean flow points.
344	! If we flow into a another grid box nbsubbas(pt)we have to make sure that its hierarchy in the
345	! basin is lower.
346	!
347	IF ( ids(pt,sbl) .GT. 0 ) THEN
348	IF ( ids(pt,sbl) .EQ. src_id .OR. ids(pt,sbl) .GT. invented_basins) THEN
349	nbfbas =nbfbas + 1
350	tbas(nbfbas) = sbl
351	tbas_hierarchy(nbfbas) = hierarchy(pt,sbl)
352	tbas_flowdir(nbfbas) = flowdir(pt,sbl)
353	ENDIF
354	ELSE
355	IF ( COUNT(coastals(pt,1:nbcoastals(pt)) .EQ. src_id) .GT. 0 ) THEN
356	nbfbas =nbfbas + 1
357	tbas(nbfbas) = sbl
358	tbas_hierarchy(nbfbas) = hierarchy(pt,sbl)
359	ENDIF
360	ENDIF
361	ENDDO
362	!
363	! Sort the possible basins by hierarchy so we can find the best placed sub-basin.
364	!
365	IF (nbfbas .GT. 1) THEN
366	!
367	IF ( MAXVAL(hierarchy(pt,1:nbsubbas(pt))) >= undef_int ) THEN
368	CALL ipslerr_p(3,'routing_bestsubbasin', &
369	"Cannot sort the sub-basins by hierarchy.", &
370	& "Some alternate method needs to be found or undef_int increased.", "")
371	ENDIF
372	!
373	DO i=1,nbfbas
374	ff = MINLOC(tbas_hierarchy(1:nbfbas))
375	fbas(i) = tbas(ff(1))
376	fbas_hierarchy(i) = tbas_hierarchy(ff(1))
377	fbas_flowdir(i) = tbas_flowdir(ff(1))
378	tbas_hierarchy(ff(1)) = undef_int
379	ENDDO
380	ELSE
381	fbas(:) = tbas(:)
382	fbas_hierarchy(:) = tbas_hierarchy(:)
383	fbas_flowdir(:) = tbas_flowdir(:)
384	ENDIF
385	!
386	! Analyse our options if there are at least one
387	!
388	sbl = undef_int
389	IF (nbfbas .EQ. 1) THEN
390	!
391	! We check that this only option is valid to provide a quality flag.
392	! If hierachy is lower then we have no problems.
393	! If the hierrachies are equal then it is possible if both sub-basins
394	! flow into the same direction.
395	!
396	IF ( fbas_hierarchy(1) < src_hierarchy ) THEN
397	sbl = fbas(1)
398	ELSE IF ( (fbas_hierarchy(1) - src_hierarchy ) < EPSILON(src_hierarchy) ) THEN
399	IF ( fbas_flowdir(1) == src_flowdir ) THEN
400	sbl = fbas(1)
401	ELSE
402	sbl = undef_int
403	ENDIF
404	ELSE
405	sbl = undef_int
406	ENDIF
407	!
408	ELSE IF (nbfbas .GT. 1) THEN
409	!
410	! First work on the hierarchy based on the sorting done above.
411	!
412	ff = MINLOC(ABS(fbas_hierarchy(1:nbfbas)-src_hierarchy))
413	!
414	! Get the basin with the hierarchy just below src_hierarchy
415	!
416	IF ( fbas_hierarchy(ff(1)) < src_hierarchy ) THEN
417	sbl = fbas(ff(1))
418	ELSE
419	!
420	!
421	IF ( ff(1) > 1 ) THEN
422	! If we are in the case of fbas_hierarchy(ff(1)) = src_hierarch we know we have a point with lower
423	! hierarchy so we flow there.
424	sbl = fbas(ff(1)-1)
425	nbopt=0
426	ELSE
427	! We are probably in the fbas_hierarchy(ff(1)) = src_hierarch situation and we need to know if we have
428	! more sub-basins in this case
429	IF ( ABS(fbas_hierarchy(ff(1))-src_hierarchy) < EPSILON(src_hierarchy) ) THEN
430	sbl = fbas(ff(1))
431	nbopt=0
432	DO i=1,nbfbas
433	IF ( ABS(fbas_hierarchy(i)-src_hierarchy) < EPSILON(src_hierarchy) ) THEN
434	nbopt = nbopt+1
435	sboptions(nbopt) = fbas(i)
436	ENDIF
437	ENDDO
438	ELSE
439	!!$ CALL ipslerr_p(2,'routing_bestsubbasin', &
440	!!$ & "The sub-basins with the closest hierrachy has still a higher value.",&
441	!!$ & "This is not an option and we will try elsewhere.", "")
442	sbl = undef_int
443	ENDIF
444	ENDIF
445
446	ENDIF
447	ELSE
448	sbl = undef_int
449	ENDIF
450	!
451	! Now we have either found a basin or we have a number of options we can explore
452	!
453	IF ( sbl == undef_int ) THEN
454	outbas = undef_int
455	quality = 0
456	ELSE
457	IF ( nbopt == 0 ) THEN
458	outbas = sbl
459	quality = 1
460	ELSE
461	!
462	! Go through the list of options and try to find a solution
463	!
464	! a) Similar flow directions. Here we are working with the integer values provided on the
465	! regular lat/lon grid. Thus we always have 8 neighbours.
466	!
467	angle(:) = undef_int
468	subbas(:) = undef_int
469	DO i=1,nbopt
470	angle(i) = MIN(MOD(src_flowdir+8 - flowdir(pt,sboptions(i)), 8), &
471	& MOD(flowdir(pt,sboptions(i))+8 - src_flowdir, 8))
472	subbas(i) = sboptions(i)
473	ENDDO
474	ff=MINLOC(angle)
475	IF (angle(ff(1)) <= 1) THEN
476	outbas = subbas(ff(1))
477	quality = 0.5
478	ENDIF
479	!
480	! b) Is the next outflow options is the ocean, then it is OK as well and even a little better.
481	!
482	DO i=1,nbopt
483	IF ( flowdir(pt,sboptions(i)) < 0 ) THEN
484	outbas = sboptions(i)
485	quality = 0.75
486	ENDIF
487	ENDDO
488	ENDIF
489	ENDIF
490	!
491	END SUBROUTINE routing_bestsubbasin
492	!
493	!
494	!!
495	!! ================================================================================================================================
496	!! SUBROUTINE : routing_reg_bestsubbasin()
497	!!
498	!>\BRIEF This routine help routing_linkup to find the most suitable sub-basin in the grid box
499	!! in which the source sub-basin flows.
500	!!
501	!! DESCRIPTION (definitions, functional, design, flags) : None
502	!!
503	!! RECENT CHANGE(S): None
504	!!
505	!! MAIN OUTPUT VARIABLE(S):
506	!!
507	!! REFERENCES : None
508	!!
509	!! FLOWCHART : None
510	!! \n
511	!_ ================================================================================================================================
512	!
513	SUBROUTINE routing_reg_bestsubbasin(src_grid, src_basin, src_id, src_hierarchy, src_fac, src_flowdir, invented_basins, &
514	& ptmax, basmax, pt, nbsubbas, ids, hierarchy, fac, flowdir, nbcoastals, coastals, &
515	& outbas, quality)
516	!
517	! Input
518	!
519	INTEGER(i_std), INTENT(in) :: src_grid, src_id, src_flowdir, src_basin
520	REAL(r_std), INTENT(in) :: src_hierarchy, src_fac ! src_fac should be integer but nevermind
521	REAL(r_std), INTENT(in) :: invented_basins
522	!
523	INTEGER(i_std), INTENT(in) :: ptmax, basmax, pt
524	INTEGER(i_std), INTENT(in) :: nbsubbas(ptmax)
525	INTEGER(i_std), INTENT(in) :: ids(ptmax,basmax), flowdir(ptmax,basmax)
526	REAL(r_std), INTENT(in) :: hierarchy(ptmax,basmax)
527	REAL(r_std), INTENT(in) :: fac(ptmax,basmax) ! fac should be integer but real is still ok
528	INTEGER(i_std), INTENT(in) :: nbcoastals(ptmax)
529	INTEGER(i_std), INTENT(in) :: coastals(ptmax,basmax)
530	!
531	! Output
532	!
533	INTEGER(i_std), INTENT(out) :: outbas
534	REAL(r_std), INTENT(out) :: quality
535	!
536	! Local
537	!
538	INTEGER(i_std) :: i, j, sbl, nbfbas, nbopt
539	INTEGER(i_std), DIMENSION(1) :: ff
540	INTEGER(i_std), DIMENSION(basmax) :: fbas, tbas, sboptions, angle, subbas
541	INTEGER(i_std), DIMENSION(basmax) :: fbas_flowdir, tbas_flowdir
542	REAL(r_std), DIMENSION(basmax) :: fbas_hierarchy, tbas_hierarchy
543	REAL(r_std), DIMENSION(basmax) :: fbas_fac, tbas_fac
544	!
545	! Debug
546	!
547	LOGICAL, PARAMETER :: debug = .FALSE. !! Logical variable to debug
548	INTEGER(i_std) :: testbasinid
549	REAL(r_std) :: maxhiediff
550	!
551	! Set some default values
552	!
553	fbas = 0
554	nbfbas = 0
555	nbopt = 0
556	sboptions(:) = undef_int
557	outbas = undef_int
558	testbasinid = 107448
559	!
560	! 1.0 Find possible basins in the list of the target grid
561	!
562	DO sbl=1,nbsubbas(pt)
563	!
564	! Either it is a standard basin or one aggregated from ocean flow points.
565	! If we flow into a another grid box nbsubbas(pt) we have to make sure that its hierarchy in the
566	! basin is lower.
567	! Trung: we get also information about flow accumulation (fac)
568	!
569	IF ( ids(pt,sbl) .GT. 0 ) THEN
570	IF ( ids(pt,sbl) .EQ. src_id .OR. ids(pt,sbl) .GT. invented_basins) THEN
571	nbfbas =nbfbas + 1
572	tbas(nbfbas) = sbl
573	tbas_hierarchy(nbfbas) = hierarchy(pt,sbl)
574	tbas_fac(nbfbas) = fac(pt,sbl)
575	tbas_flowdir(nbfbas) = flowdir(pt,sbl)
576	ENDIF
577	ELSE
578	IF ( COUNT(coastals(pt,1:nbcoastals(pt)) .EQ. src_id) .GT. 0 ) THEN
579	nbfbas =nbfbas + 1
580	tbas(nbfbas) = sbl
581	tbas_hierarchy(nbfbas) = hierarchy(pt,sbl)
582	tbas_fac(nbfbas) = fac(pt,sbl)
583	ENDIF
584	ENDIF
585	ENDDO
586	!
587	! Sort the possible basins by hierarchy so we can find the best placed sub-basin.
588	!
589	IF (nbfbas .GT. 1) THEN
590	!
591	IF ( MAXVAL(hierarchy(pt,1:nbsubbas(pt))) >= undef_int ) THEN
592	CALL ipslerr_p(3,'routing_reg_bestsubbasin', &
593	"Cannot sort the sub-basins by hierarchy.", &
594	& "Some alternate method needs to be found or undef_int increased.", "")
595	ENDIF
596	!
597	! Trung: Using "hierarchy" as more important than fac but still need using
598	! "fac" to deal with "tiny" sub-basins which can twist the stream too much
599	!
600	DO i=1,nbfbas
601	ff = MINLOC(tbas_hierarchy(1:nbfbas))
602	fbas(i) = tbas(ff(1))
603	fbas_hierarchy(i) = tbas_hierarchy(ff(1))
604	fbas_fac(i) = tbas_fac(ff(1))
605	fbas_flowdir(i) = tbas_flowdir(ff(1))
606	tbas_hierarchy(ff(1)) = undef_int
607	ENDDO
608	ELSE
609	fbas(:) = tbas(:)
610	fbas_hierarchy(:) = tbas_hierarchy(:)
611	fbas_fac(:) = tbas_fac(:)
612	fbas_flowdir(:) = tbas_flowdir(:)
613	ENDIF
614	!
615	! Debug if you want
616	!
617	IF ( debug .AND. testbasinid > 0 ) THEN
618	IF ( src_id .EQ. testbasinid ) THEN
619	WRITE (numout,*) " = = = = = = S T A R T D E B U G = = = = = = "
620	WRITE (numout,*) " Bestsubbasin TEST : ", src_id, "@", src_grid, src_basin
621	WRITE (numout,*) " nbfbas = ", nbfbas, "@", pt
622	WRITE (numout,*) " src_hierarchy = ", src_hierarchy
623	WRITE (numout,*) " src_fac = ", src_fac
624	WRITE (numout,*) " = = = = = = = = = = = = = = = = = = = = = = "
625	IF ( nbfbas .GT. 0 ) THEN
626	WRITE (numout,*) " i, fbas(i), fbas_hierarchy(i), fbas_fac(i)" !, fbas_flowdir(i)"
627	DO i=1,nbfbas
628	WRITE (numout,*) i, fbas(i), fbas_hierarchy(i), fbas_fac(i) !, fbas_flowdir(i)
629	ENDDO
630	ELSE
631	WRITE (numout,*) " We did not find any suitable sub-basin !"
632	ENDIF
633	WRITE (numout,*) " = = = = = = = E N D D E B U G = = = = = = = "
634	ENDIF
635	ENDIF
636	!
637	! Analyse our options if there are at least one
638	!
639	sbl = undef_int
640	IF (nbfbas .EQ. 1) THEN
641	!
642	! We check that this only option is valid to provide a quality flag.
643	! If hierachy is lower then we have no problems.
644	! If the hierrachies are equal then it is possible if both sub-basins
645	! flow into the same direction.
646	!
647	! Trung: If we found only one suitable sub-basin, we don't care about "fac".
648	!
649	IF ( fbas_hierarchy(1) < src_hierarchy ) THEN
650	sbl = fbas(1)
651	ELSE IF ( (fbas_hierarchy(1) - src_hierarchy ) < EPSILON(src_hierarchy) ) THEN
652	IF ( fbas_flowdir(1) == src_flowdir ) THEN
653	sbl = fbas(1)
654	ELSE
655	sbl = undef_int
656	ENDIF
657	ELSE
658	sbl = undef_int
659	ENDIF
660	!
661	ELSE IF (nbfbas .GT. 1) THEN
662	!
663	! First work on the hierarchy based on the sorting done above.
664	!
665	! Trung: Instead of using MINLOC
666	! ff = MINLOC(ABS(fbas_hierarchy(1:nbfbas)-src_hierarchy))
667	! We need a loop here: There will be situations in which
668	! there are lots of sub-basin have same ZERO hierarchy.
669	! In this case, using MINLOC will return ff(1) = 1
670	! While, in fact, it should be greater than 1.
671	!
672	ff(1) = 1
673	maxhiediff = ABS(fbas_hierarchy(ff(1))-src_hierarchy)
674	! Trung: Shorten iteration here by DO WHILE (instead of DO i = 1,nbfbas)
675	j = 0
676	i = ff(1)
677	DO WHILE ( i .GE. 1 .AND. i .LT. nbfbas .AND. j .EQ. 0 )
678	i = i + 1
679	IF ( ABS(fbas_hierarchy(i)-src_hierarchy) .LE. maxhiediff ) THEN
680	ff(1) = i
681	maxhiediff = ABS(fbas_hierarchy(ff(1))-src_hierarchy)
682	ELSE
683	j = -1
684	ENDIF
685	ENDDO
686	!
687	! Get the basin with the hierarchy just below src_hierarchy
688	!
689	! Trung: put requirement about "fac" here
690	! ONLY coastal point can have fbas_fac(ff(1)) = src_fac, right?
691	! I need to think more carefully about it.
692	! Now, I consider that there is NO fbas_fac(ff(1)) == src_fac
693	!
694	IF ( (fbas_hierarchy(ff(1)) .LT. src_hierarchy) .AND. (fbas_fac(ff(1)) .GT. src_fac) ) THEN
695	sbl = fbas(ff(1))
696	ELSE
697	!
698	!
699	IF ( ff(1) > 1 ) THEN
700	!
701	! If we are in the case of fbas_hierarchy(ff(1)) = src_hierarch we know we have a point with lower
702	! hierarchy so we flow there.
703	!
704	! Trung: We consider "fbas_fac(ff(1)) > src_fac", so we need to go
705	! higher in sorted list fbas(:) to find suitable sub-basin.
706	!
707	j = 0
708	i = ff(1)
709	DO WHILE ( i .GT. 1 .AND. j .EQ. 0 )
710	i = i - 1
711	IF ( fbas_hierarchy(i) .LT. src_hierarchy .AND. fbas_fac(i) .GT. src_fac ) THEN
712	sbl = fbas(i)
713	nbopt = 0
714	j = -1
715	ENDIF
716	ENDDO
717	ELSE
718	! We are probably in the fbas_hierarchy(ff(1)) = src_hierarch situation and we need to know if we have
719	! more sub-basins in this case
720	! Trung: If the quality of routing.nc is good, we will never come to
721	! this step, right?
722	IF ( ABS(fbas_hierarchy(ff(1))-src_hierarchy) < EPSILON(src_hierarchy) ) THEN
723	sbl = fbas(ff(1))
724	nbopt=0
725	DO i=1,nbfbas
726	! Trung: the condition about "fac" here need to be investigated later
727	IF ( (ABS(fbas_hierarchy(i)-src_hierarchy) < EPSILON(src_hierarchy)) .AND. (fbas_fac(i) > src_fac) ) THEN
728	nbopt = nbopt+1
729	sboptions(nbopt) = fbas(i)
730	ENDIF
731	ENDDO
732	ELSE
733	!!$ CALL ipslerr_p(2,'routing_reg_bestsubbasin', &
734	!!$ & "The sub-basins with the closest hierrachy has still a higher value.",&
735	!!$ & "This is not an option and we will try elsewhere.", "")
736	sbl = undef_int
737	ENDIF
738	ENDIF
739	ENDIF
740	ENDIF
741	!
742	! Now we have either found a basin or we have a number of options we can explore
743	!
744	IF ( sbl == undef_int ) THEN
745	outbas = undef_int
746	quality = 0
747	ELSE
748	IF ( nbopt == 0 ) THEN
749	outbas = sbl
750	quality = 1
751	ELSE
752	!
753	! Go through the list of options and try to find a solution
754	!
755	! a) Similar flow directions. Here we are working with the integer values provided on the
756	! regular lat/lon grid. Thus we always have 8 neighbours.
757	!
758	angle(:) = undef_int
759	subbas(:) = undef_int
760	DO i=1,nbopt
761	angle(i) = MIN(MOD(src_flowdir+8 - flowdir(pt,sboptions(i)), 8), &
762	& MOD(flowdir(pt,sboptions(i))+8 - src_flowdir, 8))
763	subbas(i) = sboptions(i)
764	ENDDO
765	ff=MINLOC(angle)
766	IF (angle(ff(1)) <= 1) THEN
767	outbas = subbas(ff(1))
768	quality = 0.5
769	ENDIF
770	!
771	! b) Is the next outflow options is the ocean, then it is OK as well and even a little better.
772	!
773	DO i=1,nbopt
774	IF ( flowdir(pt,sboptions(i)) < 0 ) THEN
775	outbas = sboptions(i)
776	quality = 0.75
777	ENDIF
778	ENDDO
779	ENDIF
780	ENDIF
781	!
782	END SUBROUTINE routing_reg_bestsubbasin
783	!! ================================================================================================================================
784	!! SUBROUTINE : routing_updateflow
785	!!
786	!>\BRIEF This subroutine allows to update the various variables which store the flow
787	!! from one gridbox to another. It is just a tool to simplify the main code.
788	!!
789	!! DESCRIPTION (definitions, functional, design, flags) : None
790	!!
791	!! RECENT CHANGE(S): None
792	!!
793	!! MAIN OUTPUT VARIABLE(S):
794	!!
795	!! REFERENCES : None
796	!!
797	!! FLOWCHART : None
798	!! \n
799	!_ ================================================================================================================================
800	!
801	SUBROUTINE routing_updateflow(sp, sb, ep, eb, pt, bas, basmax, out_grid, out_basin, in_number, in_grid, in_basin)
802	!
803	! ARGUMENTS
804	!
805	INTEGER(i_std), INTENT(in) :: sp, sb !! Source point and source basin
806	INTEGER(i_std), INTENT(in) :: ep, eb !! target point and target basin
807	INTEGER(i_std), INTENT(in) :: pt, bas, basmax !! dimensions of arrays
808	!
809	INTEGER(i_std), DIMENSION(pt,bas), INTENT(inout) :: out_grid !! Grid to which we flow
810	INTEGER(i_std), DIMENSION(pt,bas), INTENT(inout) :: out_basin !! Basin in the out_grid
811	INTEGER(i_std), DIMENSION(pt,basmax), INTENT(inout) :: in_number !! Number of basins flowing into this grid
812	INTEGER(i_std), DIMENSION(pt,basmax,basmax), INTENT(inout) :: in_grid !! Source grid for this grid
813	INTEGER(i_std), DIMENSION(pt,basmax,basmax), INTENT(inout) :: in_basin !! Source basin for this grid
814	!
815	! LOCAL
816	!
817	INTEGER(i_std) :: i
818	LOGICAL :: transfer
819	!
820	! Test that the target sub-basin is not within the sources (inflows) of the current point
821	!
822	transfer = .TRUE.
823	DO i=1,in_number(sp,sb)
824	IF (in_grid(sp,sb,i) == ep .AND. in_basin(sp,sb,i) == eb ) THEN
825	WRITE(numout,*) "We have a problem here !", ep, eb, " already flows into ", sp, sb
826	transfer = .FALSE.
827	ENDIF
828	ENDDO
829	!
830	! Just in case this was not caught in the calling program.
831	!
832	IF ( ep <= 0 .OR. ep == undef_int .OR. eb <= 0 .OR. eb == undef_int) THEN
833	transfer = .FALSE.
834	ENDIF
835	!
836	IF ( transfer ) THEN
837	out_grid(sp,sb) = ep
838	out_basin(sp,sb) = eb
839	!
840	in_number(ep,eb) = in_number(ep,eb) + 1
841	IF ( in_number(ep,eb) .LE. basmax ) THEN
842	in_grid(ep,eb,in_number(ep,eb)) = sp
843	in_basin(ep,eb,in_number(ep,eb)) = sb
844	ELSE
845	WRITE(numout,*) "routing_updateflow : Too many basins flowing into this point ", ep, eb
846	CALL ipslerr_p(3,'routing_updateflow', &
847	"Too many basins flowing into this point ",'Increase nbxmax or its derived quantities','')
848	ENDIF
849	ENDIF
850	!
851	END SUBROUTINE routing_updateflow
852	!
853	!! ================================================================================================================================
854	!! SUBROUTINE : routing_sortcoord
855	!!
856	!>\BRIEF This subroutines orders the coordinates to go from North to South and West to East.
857	!!
858	!! DESCRIPTION (definitions, functional, design, flags) : None
859	!!
860	!! RECENT CHANGE(S): None
861	!!
862	!! MAIN OUTPUT VARIABLE(S):
863	!!
864	!! REFERENCES : None
865	!!
866	!! FLOWCHART : None
867	!! \n
868	!_ ================================================================================================================================
869
870	SUBROUTINE routing_sortcoord(nb_in, coords, direction, nb_out)
871	!
872	IMPLICIT NONE
873	!
874	!! INPUT VARIABLES
875	INTEGER(i_std), INTENT(in) :: nb_in !!
876	REAL(r_std), INTENT(inout) :: coords(nb_in) !!
877	!
878	!! OUTPUT VARIABLES
879	INTEGER(i_std), INTENT(out) :: nb_out !!
880	!
881	!! LOCAL VARIABLES
882	CHARACTER(LEN=2) :: direction !!
883	INTEGER(i_std) :: ipos !!
884	REAL(r_std) :: coords_tmp(nb_in) !!
885	INTEGER(i_std), DIMENSION(1) :: ll !!
886	INTEGER(i_std) :: ind(nb_in) !!
887
888	!_ ================================================================================================================================
889	!
890	ipos = 1
891	nb_out = nb_in
892	!
893	! Compress the coordinates array
894	!
895	DO WHILE ( ipos < nb_in )
896	IF ( coords(ipos+1) /= undef_sechiba) THEN
897	IF ( COUNT(coords(ipos:nb_out) == coords(ipos)) > 1 ) THEN
898	coords(ipos:nb_out-1) = coords(ipos+1:nb_out)
899	coords(nb_out:nb_in) = undef_sechiba
900	nb_out = nb_out - 1
901	ELSE
902	ipos = ipos + 1
903	ENDIF
904	ELSE
905	EXIT
906	ENDIF
907	ENDDO
908	!
909	! Sort it now
910	!
911	! First we get ready and adjust for the periodicity in longitude
912	!
913	coords_tmp(:) = undef_sechiba
914	IF ( INDEX(direction, 'WE') == 1 .OR. INDEX(direction, 'EW') == 1) THEN
915	IF ( MAXVAL(ABS(coords(1:nb_out))) .GT. 160 ) THEN
916	coords_tmp(1:nb_out) = MOD(coords(1:nb_out) + 360.0, 360.0)
917	ELSE
918	coords_tmp(1:nb_out) = coords(1:nb_out)
919	ENDIF
920	ELSE IF ( INDEX(direction, 'NS') == 1 .OR. INDEX(direction, 'SN') == 1) THEN
921	coords_tmp(1:nb_out) = coords(1:nb_out)
922	ELSE
923	WRITE(numout,*) 'The chosen direction (', direction,') is not recognized'
924	CALL ipslerr_p(3,'routing_sortcoord','The chosen direction is not recognized','First section','')
925	ENDIF
926	!
927	! Get it sorted out now
928	!
929	ipos = 1
930	!
931	IF ( INDEX(direction, 'WE') == 1 .OR. INDEX(direction, 'SN') == 1) THEN
932	DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1)
933	ll = MINLOC(coords_tmp(:), coords_tmp /= undef_sechiba)
934	ind(ipos) = ll(1)
935	coords_tmp(ll(1)) = undef_sechiba
936	ipos = ipos + 1
937	ENDDO
938	ELSE IF ( INDEX(direction, 'EW') == 1 .OR. INDEX(direction, 'NS') == 1) THEN
939	DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1)
940	ll = MAXLOC(coords_tmp(:), coords_tmp /= undef_sechiba)
941	ind(ipos) = ll(1)
942	coords_tmp(ll(1)) = undef_sechiba
943	ipos = ipos + 1
944	ENDDO
945	ELSE
946	WRITE(numout,*) 'The chosen direction (', direction,') is not recognized (second)'
947	CALL ipslerr_p(3,'routing_sortcoord','The chosen direction is not recognized','Second section','')
948	ENDIF
949	!
950	coords(1:nb_out) = coords(ind(1:nb_out))
951	IF (nb_out < nb_in) THEN
952	coords(nb_out+1:nb_in) = zero
953	ENDIF
954	!
955	END SUBROUTINE routing_sortcoord
956	!
957	!! ================================================================================================================================
958	!! SUBROUTINE : routing_diagbox
959	!!
960	!>\BRIEF This function will determine if this point is to be diagnosed or not.
961	!!
962	!! DESCRIPTION (definitions, functional, design, flags) : None
963	!!
964	!! RECENT CHANGE(S): None
965	!!
966	!! MAIN OUTPUT VARIABLE(S):
967	!!
968	!! REFERENCES : None
969	!!
970	!! FLOWCHART : None
971	!! \n
972	!_ ================================================================================================================================
973	!
974	LOGICAL FUNCTION routing_diagbox(ip, diaglalo)
975	!
976	INTEGER(i_std), INTENT(in) :: ip
977	REAL(r_std), DIMENSION(:,:), INTENT(in) :: diaglalo
978	!
979	! Local
980	!
981	REAL(r_std) :: dlon, dlat
982	INTEGER(i_std) :: ij, nbdiag
983	!
984	nbdiag=SIZE(diaglalo,DIM=1)
985	!
986	dlon=ABS(corners(ip,NbSegments,1)-corners(ip,1,1))
987	dlat=ABS(corners(ip,NbSegments,2)-corners(ip,1,2))
988	DO ij=1,NbSegments-1
989	dlon = MAX(dlon,ABS(corners(ip,ij+1,1)-corners(ip,ij,1)))
990	dlat = MAX(dlat,ABS(corners(ip,ij+1,2)-corners(ip,ij,2)))
991	ENDDO
992	!
993	routing_diagbox = .FALSE.
994	DO ij=1,nbdiag
995	IF ( ABS(lalo(ip,1)-diaglalo(ij,1)) < dlat/2.0 .AND. &
996	& ABS(lalo(ip,2)-diaglalo(ij,2)) < dlon/2.0 ) THEN
997	routing_diagbox = .TRUE.
998	ENDIF
999	ENDDO
1000	!
1001	END FUNCTION routing_diagbox
1002	!
1003	LOGICAL FUNCTION routing_diagbox_g(ip, diaglalo)
1004	!
1005	INTEGER(i_std), INTENT(in) :: ip
1006	REAL(r_std), DIMENSION(:,:), INTENT(in) :: diaglalo
1007	!
1008	! Local
1009	!
1010	REAL(r_std) :: dlon, dlat
1011	INTEGER(i_std) :: ij, nbdiag
1012	!
1013	IF (.NOT. is_root_prc) THEN
1014	WRITE(,) "The function routing_diagbox_g is written to work on the root proc and the full grid"
1015	STOP "routing_diagbox_g"
1016	ENDIF
1017	!
1018	nbdiag=SIZE(diaglalo,DIM=1)
1019	!
1020	dlon=ABS(corners_g(ip,NbSegments,1)-corners_g(ip,1,1))
1021	dlat=ABS(corners_g(ip,NbSegments,2)-corners_g(ip,1,2))
1022	DO ij=1,NbSegments-1
1023	dlon = MAX(dlon,ABS(corners_g(ip,ij+1,1)-corners_g(ip,ij,1)))
1024	dlat = MAX(dlat,ABS(corners_g(ip,ij+1,2)-corners_g(ip,ij,2)))
1025	ENDDO
1026	!
1027	routing_diagbox_g = .FALSE.
1028	DO ij=1,nbdiag
1029	IF ( ABS(lalo_g(ip,1)-diaglalo(ij,1)) < dlat/2.0 .AND. &
1030	& ABS(lalo_g(ip,2)-diaglalo(ij,2)) < dlon/2.0 ) THEN
1031	routing_diagbox_g = .TRUE.
1032	ENDIF
1033	ENDDO
1034	!
1035	END FUNCTION routing_diagbox_g
1036	!
1037	END MODULE routing_tools

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

Download in other formats: